Minh Son Phan

Yukiko Kenmochi

Akihiro Sugimoto

Hugues Talbot

Eric Andres

Rita Zrour

Efficient robust digital annulus fitting with bounded error

Keywords: fitting, annulus, approximation, halfspace range searching

A digital annulus is defined as a set of grid points lying between two circles sharing an identical center and separated by a given width. This paper deals with the problem of fitting a digital annulus to a given set of points in a 2D bounded grid. More precisely, we tackle the problem of finding a digital annulus that contains the largest number of inliers. As the current best algorithm for exact optimal fitting has a computational complexity in O(N 3 log N ) where N is the number of grid points, we present an approximation method featuring linear time complexity and bounded error in annulus width, by extending the approximation method previously proposed for digital hyperplane fitting. Experiments show some results and runtime in practice.

Introduction

Shape fitting is widely used in computer vision, computational geometry, image analysis and many other areas. Given a set of points S, a shape model M (for example, lines, curves, planes, or circles), and some criterion F , the problem is to fit M to S under F . Depending on the criteria employed, we have several well-known methods dealing with this problem. The method of least-squares [START_REF] Boyd | Convex optimization[END_REF] minimizes the sum of squared residuals from all given points. This method is popular, but it is not efficient in the presence of outliers. In order to enhance the robustness, many alternative methods have also been proposed, such as the method of least-absolute-values [START_REF] Press | Numerical recipes: the art of scientific computing[END_REF], the method of least-median-of-squares [START_REF] Rousseeuw | Least median of squares regression[END_REF]. However, these methods are still not robust in the presence of many outliers. One of the commonly used methods for robust shape fitting is RANdom SAmple Consensus (RANSAC) [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF]. It can provide robust results even if there is a significant number of outliers, with reasonable runtime in practice, while it provides neither guarantee of optimal solution nor bound of approximation error.

In this paper, we use a criterion similar to that of RANSAC. As we assume that our input data is discrete, such as a digital image, we use a digital shape model, instead of a continuous one as generally used for the above mentioned conventional methods. Fitting of digital lines [START_REF] Zrour | Optimal consensus set for digital line and plane fitting[END_REF] or digital planes [START_REF] Aiger | Efficient robust digital hyperplane fitting with bounded error[END_REF] has been already studied. This paper aims at exploring the problem of fitting a digital circle. We consider a digital circle as a digital annulus [START_REF] Andres | The Discrete Analytical Hyperspheres[END_REF], defined as a set of grid points lying between two circles sharing a center but with different radii. This paper treats the problem of fitting a digital annulus to a given set of points in a 2D bounded grid. More precisely, we tackle the problem of finding a digital annulus that contains the largest number of inliers. The current best algorithm for exact optimal fitting has a computational complexity in O(N 3 log N ) [START_REF] Largeteau-Skapin | O(n 3 log n) time complexity for the optimal consensus set computation for Andres and 4-connected Flake circles[END_REF]. We present an approximation method featuring linear time complexity and bounded error in annulus width by extending an approximation method proposed for digital hyperplane fitting [START_REF] Aiger | Efficient robust digital hyperplane fitting with bounded error[END_REF].

We first formulate our digital annulus fitting as the digital plane fitting problem with some particular setting. For our approximation method, we then employ the approximate halfspace range searching technique [START_REF] Da Fonseca | Approximate range searching: The absolute model[END_REF], which consists of two main steps: generating an approximate range space, and providing a data structure to solve approximate halfspace range searching with a constant query time. We follow these two steps, in the context of solving our converted digital plane fitting, by combining two problems of approximate halfspace range searching. We also explain how to interpret results of digital plane fitting in the context of digital annulus fitting, with an interpretation of bounded error. Finally, we provide experimental results and measured runtimes.

Problem of digital annulus fitting

Digital annulus and its fitting problem

In this paper, we follow the definition of a digital annulus as introduced by Andres et al. [START_REF] Andres | The Discrete Analytical Hyperspheres[END_REF] and used by Zrour et al. for fitting [START_REF] Largeteau-Skapin | O(n 3 log n) time complexity for the optimal consensus set computation for Andres and 4-connected Flake circles[END_REF]. We first describe the continuous model of circles. A continuous circle C with center (a, b) and radius r in the Euclidean space R 2 is defined by

C = {(x, y) ∈ R 2 : (x -a) 2 + (y -b) 2 = r 2 } (1)
where a, b, r ∈ R. In the discrete space Z 2 , the digitization of ( 1) is defined by a set of grid points lying between two circles that share the same center (a, b) and two different radii r -w 2 and r + w 2 , where w ∈ R is the width between the two circles. Namely,

A = (x, y) ∈ Z 2 : r - w 2 2 ≤ (x -a) 2 + (y -b) 2 ≤ r + w 2 2 . ( 2
)
This is called a digital annulus, or Andres circle [START_REF] Andres | The Discrete Analytical Hyperspheres[END_REF]. Given a set S of discrete points coming from a [0, δ] 2 grid, where δ ∈ N, the problem is to fit this digital annulus model A to S such that the fitted digital annulus contains as many points of S as possible, called inliers. This problem is described below. Points that are not contained in the model are called outliers.

Problem 1 (Digital annulus fitting) Input : A set S of discrete points bounded in a [0, δ] 2 grid, and a fixed width w.

Output : A digital annulus A containing the maximum number of inliers in S.

This digital annulus fitting problem is solved by an exact method, introduced in [START_REF] Largeteau-Skapin | O(n 3 log n) time complexity for the optimal consensus set computation for Andres and 4-connected Flake circles[END_REF], in time O(|S| 3 log |S|). This complexity of the exact solution is not feasible in practice when |S| is large. In some cases it may be sufficient to compute an approximate solution if the solution is obtained with inexpensive computational complexity. We thus propose here an approximate method with linear runtime and bounded error, instead of searching for the exact solution. For this purpose, we first convert our digital annulus fitting problem into a more simple problem in the following.

Digital plane fitting induced by digital annulus fitting

Since the inequalities (2) are non-linear with respect to the unknown parameters a, b, r, we convert those to linear ones along the known conversion [START_REF] O'rourke | Computing circular separability[END_REF]. Then, the digital annulus formula of (2) can be converted into the digital plane formula [START_REF] Andres | Discrete Analytical Hyperplanes[END_REF], defined by two inequalities of the form 0 ≤ Ax + By + z + C ≤ W where W is the width.

The inequalities (2) are written as:

r - w 2 2 ≤ -2ax -2by + x 2 + y 2 + a 2 + b 2 ≤ r + w 2 2 . ( 3
) Letting            z = x 2 + y 2 , A = -2a, B = -2b, C = a 2 + b 2 -r -w 2 2 , C = a 2 + b 2 -r + w 2 2 , (4)
we have that 0

≤ Ax + By + z + C ≤ C -C .
We can thus define the converted digital plane model as follows:

P = (x, y, z) ∈ Z 3 : 0 ≤ Ax + By + z + C ≤ C -C . (5)
Obviously this has the form of a digital plane [START_REF] Andres | Discrete Analytical Hyperplanes[END_REF], which is defined as a set of grid points lying between two parallel planes, whose width is set to be W = C -C . Therefore, (5) leads to the following observation.

Observation 1 The problem of digital annulus fitting is equivalent to the problem of digital plane fitting for the model defined by (5). This converted problem is described as follows.

Problem 2 (Converted digital plane fitting) Input : A set S of discrete points bounded in a [0, δ] 2 grid, and width w.

Output : A digital plane P containing the maximum number of inliers when the input set is S = {(x, y, x 2 + y 2 ) : (x, y) ∈ S}.

Figure 1 illustrates the relationship between A and P, so that the 2D coordinates (x, y) of A are converted into 3D coordinates (x, y, x 2 + y 2 ) of P. In fact, the setting of new parameters in (4) is not efficient since ∀(x, y) ∈ [0, δ] 2 ⇒ z = x 2 + y 2 ∈ [0, 2δ 2 ], and consequently the 2D grid of size [0, δ] 2 for digital annulus fitting is changed to the 3D grid of size [0, δ] 2 × [0, 2δ 2 ] for the converted digital plane fitting; the grid size for z is too large in practice. We therefore define, instead of (4),

               z = 1 2δ (x 2 + y 2 ), A = -a δ , B = -b δ , C = 1 2δ a 2 + b 2 -r -w 2 2 , C = 1 2δ a 2 + b 2 -r + w 2 2 , (6)
by dividing all the terms in (3) by 2δ. With this setting, we see that

x, y ∈ [0, δ] ⇒ z = 1 2δ (x 2 + y 2 ) ∈ [0, δ]
, so that we have a 3D grid [0, δ] 3 for digital plane fitting.

Concerning the ranges of the new parameters A, B, C and C , we first obtain

a, b ∈ [0, δ] ⇒ A = -a δ , B = -b δ ∈ [-1, 0]
. Accordingly, the search space for the converted digital plane is reduced, compared with the standard digital plane fitting. To enforce r -w 2 2 and r + w 2 2 to be positive, C and C in (6) must satisfy

C, C ≤ a 2 + b 2 2δ or equivalently, C, C ≤ δ(A 2 + B 2 ) 2 .
This indicates that the ranges of the parameters C and C depend on those of A and B, and we can simply specify that

C, C ∈ [-δ, δ(A 2 +B 2 ) 2
]. We also recall that the width W of a digital plane, converted from a digital annulus, is given by

W = C -C = rw δ (7)
after the normalization of ( 6). This indicates that W is not constant, but depends on the radius r of the digital annulus. Therefore, our converted digital plane fitting problem is different from the standard one [START_REF] Aiger | Efficient robust digital hyperplane fitting with bounded error[END_REF] in that the width W is not a constant. Since we obtain from ( 6)

r = δ(A 2 + B 2 ) 2 -2δC + w 2 , (8)
we can also say that W depends on A, B and C. In other words, once we set the values for A, B and C, we can automatically set the value for r, and thus for W as well.

For this converted digital plane fitting problem, we apply the approximate method, which is originally proposed for the standard digital plane fitting problem [START_REF] Aiger | Efficient robust digital hyperplane fitting with bounded error[END_REF]. We show in the next section that the dependency of W on A, B and C does not critically obstruct applying the approximate method to our converted digital plane fitting problem.

3 Approximate method for digital annulus fitting

Method outline

The proposed digital annulus fitting algorithm consists of four steps, as described in Algorithm 1. In the first step, we convert the problem of digital annulus fitting into a problem of digital plane fitting by using the notions presented in the previous section. In other words, a set of points S in the 2D grid of size [0, δ] 2 is transformed into the set of points S = {(x, y, x 2 +y 2 2δ ) : (x, y) ∈ S} in the 3D grid of size [0, δ] 3 . In the second step, we build a data structure with respect to S for the approximate halfspace range searching algorithm. More precisely, we specify a query set H for approximate halfspace range searching with a given bounded error ε, and generate a data structure that allows to answer the following query: for any plane h ∈ H, how many points of S are on or below h with the bounded error ε. In the third step, we first modify the query for the digital plane fitting with width W + 5ε, i.e., we count the number of points of S that lie between two parallel planes with distance W + 5ε. We then find an approximate digital plane with the maximum number of inliers. The final step is to interpret the result of the converted problem as that of the original digital annulus problem.

As Steps 1 and 4 are already explained in the previous section, we describe Steps 2 and 3. Furthermore, we discuss the complexity of our proposed method.

Algorithm 1: Approximate digital annulus fitting

Input: A set S of discrete points bounded in [0, δ] 2 grid, an approximation error ε > 0, a digital annulus width w > 0. Output: The fitted digital annulus A containing the maximum number of inliers in S. 1 Convert Problem 1 into Problem 2; 2 Build a data structure for approximate halfspace range searching; 3 Make queries about range counting for all of the approximate digital planes and Find the one that contains the maximum number of inliers; 4 Interpret the result of the digital plane fitting as digital annulus fitting.

Data structure for approximate halfspace range searching

As this part is based on the work of Fonseca et al. [START_REF] Da Fonseca | Approximate range searching: The absolute model[END_REF], we first summarize their method, and then explain the changes suitable for our problem.

Given a set of 3D points T and a plane h in the form Ax + By + z + C = 0, halfspace range searching (or counting) is the problem of counting the number of points of T that are on, or below h. In order to solve this problem, Fonseca et al. [START_REF] Da Fonseca | Approximate range searching: The absolute model[END_REF] have proposed an approximate method with a bounded error ε > 0.

One of the main ideas of the method is to define a sufficiently large finite set of planes H, so that any query plane is approximated by some plane in H with bounded error ε. In this paper, H is obtained by generating different slopes

A, B ∈ [-1, 0] and different intercepts C ∈ [-δ, δ(A 2 +B 2 )

2

]; for their ranges, see the previous section. For A and B, the interval of their variations should be set to be ε δ , while it should be ε for C (see [START_REF] Da Fonseca | Approximate range searching: The absolute model[END_REF] for more details). For each plane h ∈ H, we compute a counting function f (T ∩ h) that returns approximately the number of points of T that are on, or below h. The objective of constructing a data structure is to compute this function f (T ∩ h) for all the planes h ∈ H. Points within distance ε from a plane h may or may not be reported.

The computational complexity is described as follows:

Theorem 1. (Fonseca et al. [START_REF] Da Fonseca | Approximate range searching: The absolute model[END_REF]) For a set of N points in the [0, 1] 3 unit grid and some ε > 0, one can build a data structure with O(ε -3 ) storage space, in O(N + ε -3 log O(1) (ε -1 )) time, such that for a given query plane h, the number of points on or below h can be approximately reported in O(1) time, in the following sense: all the points (below h) that have a larger distance than ε from h are counted. Points that are closer to h on both sides may or may not be reported.

In fact, in order to obtain this complexity (linear to N ), we build a data structure based on an octree. Let us call a [0, δ] 3 grid a primitive cube, which is divided into 8 children. Its children are also divided recursively into 8 subcubes, until the size of divided cubes equals ε. Let X be a cube generated from a primitive cube, and T be a set of points bounded by X in dimension δ. We then verify whether each plane of H passes through X, and compute f (T ∩ h), ∀h ∈ H in the following way. Let X i denote the children of X, for i = 1..8, and T i = X i ∩ T . We then compute f (T ∩ h) by summarizing all the results from its children f (T i ∩ h), which is also recursively computed, as follows:

f (T ∩ h) =        8 i=1 f (T i ∩ h) if the size of X is larger than ε and h ∩ X = ∅, 0
if the size of X is larger than ε and h is below X, |T | otherwise. ( 9) Finally, the data structure of a primitive cube is built recursively. Figure 2 illustrates a 2D example for computing the counting function f (T ∩ h) of ( 9).

Query for approximate digital plane fitting

After building a data structure for a set of points S in the [0, δ] 3 grid, we count the number of points lying in every digital plane, i.e., between a pair of parallel planes with distance W . In order to calculate it approximately, we need a finite query set Q of approximate digital planes for our converted digital plane fitting as follows. For each plane h ∈ H, we take the values of A, B and C of h for the digital plane parameters A, B, and C, respectively. Once the values of A, B and C are fixed, we can automatically set C or W from ( 7) and ( 8), i.e., the rest of the digital plane parameters. To generate a query for approximate digital plane fitting with a width W , we refer to the following theorem presented in [START_REF] Aiger | Efficient robust digital hyperplane fitting with bounded error[END_REF]. 5Theorem 2. (Aiger et al. [START_REF] Aiger | Efficient robust digital hyperplane fitting with bounded error[END_REF]) Given a set of N points in [0, δ] 3 , and some ε > 0, W > 0, a digital plane of width W + 5ε that contains n > n opt points, can be found in O(N + ( δ ε ) 3 log O(1) ( δ ε )) time, where n opt is the maximum number of points that any digital plane of width W in [0, δ] 3 can contain.

In fact, we need to generate a new query set Q such that any digital plane with width W is completely contained in at least one of the digital planes in Q. Theorem 2 indicates that such query digital planes must have a width of at least W + 5ε. In other words, due to this setting, the important property, n ≥ n opt , is guaranteed. Even if W is not constant as in this paper, because it depends on each parameter set of A, B and C, this result is still valid.

An approximate solution of the converted digital plane fitting is obtained by finding the digital plane in Q in which the number of points is maximized. 3. Relationship between approximate (dotted line) and optimal digital line (continuous line). There is at least one approximate digital line that contains the optimal digital line.

Complexity of the approximate digital annulus fitting

As a result, we obtain the next corollary from Theorem 2 after minor adaptations.

Corollary 1 Given a set of N points on a grid [0, δ] 2 , and some ε > 0, w > 0, a digital annulus of width w + 5 δ r ε that contains n > n opt points, can be found in

O(N + ( δ ε ) 3 log O(1) ( δ ε ))
time, where r is the radius of a digital annulus of width w + 5 δ r ε and n opt is the maximum number of points that any digital annulus of width w in [0, δ] 2 can contain.

This corollary shows that the proposed method has linear computational complexity with respect to N . Obviously, if we decrease the bounded error ε, the computational time increases due to the second term, just like for approximated digital plane fitting; see Theorem 2. The difference to Theorem 2 is that ε cannot be tuned directly in the original 2D grid for digital annulus fitting, but in the 3D grid for the converted digital plane fitting. In other words, we can give ε as a constant ambiguous zone around the border of a digital plane in 3D (i.e. a zone in which points may or may not be counted), but once we project it to 2D before the conversion, the corresponding ambiguous zone becomes δ r ε which is obtained from [START_REF] Andres | The Discrete Analytical Hyperspheres[END_REF].

Experiments

The 2D point cloud data used for the experiments were created as follows. By using the digital annulus A of (2) with w = 3.0, a = b = 100.0, we randomly generated 400 inliers that satisfy A, and 100n outliers that do not satisfy A, for n = 0, . . . , 10; there are 11 variations for the number of outliers. All of the generated points are bounded in the [0, 200] 2 grid. The values for the parameters of the converted digital planes P of (5) are also computed: A = B = -0.5, C = 30.1975 and C = 28.84. Note that the experiments were run on a standard PC with core i3 CPU at 2.20 GHz.

In the first experiment, we had n = 1 constant (i.e., 100 outliers) for the input data and observed the fitting results by using different bounded errors such that ε = 3.0, 2.5, 2.0, 1.5, 1.0. As seen in Table 1 and Figure 4, the number of inliers decreases as ε becomes smaller, and it tends to converge to the ground-truth solution. Clearly, the smaller the value of ε, the more precise the solution. Concerning the runtime, it is in fact polynomial in factor ε (see Figure ). Therefore, we need to select an appropriate ε (for example 1.5 or 1.0 for this experiment) for assuring a reasonable runtime and an approximate solution which is not far from the ground-truth.

In the second experiment, we decided to use ε = 1.5 for digital annulus fitting, for all the input data sets (i.e., different numbers of outliers, 100n for n = 0, . . . , 10). As seen in Table 2 and Figure 5, numbers of inliers are always greater than 400 (i.e., the number of inliers for the ground-truth digital annulus) and the parameter values are almost similar, while the runtimes are linearalmost constant (around 1 minute). For comparison, the exact algorithm [START_REF] Largeteau-Skapin | O(n 3 log n) time complexity for the optimal consensus set computation for Andres and 4-connected Flake circles[END_REF] with computational complexity in O(N 3 log N ) was also applied. The obtained parameter values are a = 99.25, b = 99.5 and r = 90.65 for all n. The runtimes, shown in Table 2 and Figure 6, are indeed polynomial with respect to the number of points.

Table 2. Results of digital annulus fitting with ε = 1.5 to the data with 100n outliers. The runtimes of the approximation method and the exact method are compared.

Conclusion

This paper presented an approximation method for fitting a digital annulus.

Given N points in a [0, δ] 2 grid, a width w and an approximation parameter ε, the fitted digital annulus that contains the maximum number of points with a width w + 5 δ r ε, is found in O(N + ε -3 log O(1) (ε -1 )) time. The method is linear in the number of points N , but it is polynomial in approximation parameter ε. Therefore, we need to find an appropriate ε for assuring a reasonable runtime. The method is also robust to outliers. The number of inliers converges to the optimal solution when ε decreases.

Fig. 1 .

 1 Fig. 1. Illustration of converting a circle (colored in blue) in 2D into a plane (colored in red) in 3D. The converted 3D points are on the parabola z = x 2 + y 2 (colored in violet).

Fig. 2 .

 2 Fig. 2. Illustration for recursively computing f (T ∩ h) in a 2D case (a line) by dividing a primitive square (corresponding to a cube in 3D) into 4 children. The light grey cells are the squares of size ε through which a query line h goes, while the dark grey cells are the squares below h. See (9) for the definition.

Fig.

  Fig. 3. Relationship between approximate (dotted line) and optimal digital line (continuous line). There is at least one approximate digital line that contains the optimal digital line.

Fig. 4 .

 4 Fig. 4. The fitted digital annuli (colored in red) for ε = 1.5 (a) and ε = 1.0 (b). They contain the true digital annulus (colored in blue). The converted digital planes are also illustrated for ε = 1.5 (c) and ε = 1.0 (d).

Table

1 .

 1 Results of digital annulus fitting with varied ε to the data with n = 1. 80 sec 96.0 96.0 89.62 -0.48 -0.48 34.0 18.0 441 2.5 13.11 sec 92.5 97.5 90.89 -0.46 -0.49 31.5 17.5 435 2.0 26.43 sec 96.0 104.0 95.03 -0.48 -0.52 33.0 22.0 425 1.5 59.46 sec 99.0 96.0 90.64 -0.48 -0.50 31.5 22.5 419 1.0 200.54 sec 100.0 98.0 89.46 -0.50 -0.49 32.0 26.0 406 ground truth 100.0 100.0 90.5 -0.50 -0.50 30.20 28.84 400 6

Fig. 5 .

 5 Fig. 5. The fitted digital annuli (colored in red) with 100n outliers in input for n = 1 (a), n = 4 (b), n = 10 (c). They contains the true digital annulus (colored in blue). The converted digital planes are also illustrated for n = 1 (d), n = 4 (e), n = 10 (f).

Even if the original theorem was proposed for N points on a [0, δ] 3 grid, the theorem itself, as well as the proposed method, is established for a more general setting such that the N points can simply be in [0, δ] 3 regardless of any grid.

This work has been partly supported by the French Agence Nationale de la Recherche (ANR-2010-BLAN-0205 03) and a French-Japanese joint research project called the Sakura program (No. 27608XJ).


Julien Toulouse
email: julien.toulouse@upmc.fr

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction to the calculation of molecular properties by response theory

Introduction to the calculation of molecular properties by response theory

This document reviews the calculation of time-independent and time-dependent molecular properties in quantum chemistry. Sections 1.1, 1.2 and 1.3 are inspired to a large extent by the lectures of Prof. Trygve Helgaker (University of Olso, Norway) [START_REF] Helgaker | Time-Independent Molecular Properties[END_REF]. For more intensive reviews, see Refs. [START_REF] Mcweeny | Methods of Molecular Quantum Mechanics[END_REF][START_REF] Helgaker | [END_REF][START_REF] Saue | Relativistic Electronic Structure Theory[END_REF].

1 Time-independent molecular properties

Definition

We consider an isolated molecular system of Hamiltonian Ĥ0 , and apply a static perturbation V (x), depending on a variable x (representing the strength of the perturbation), such that the perturbation vanishes for x = 0. The perturbed Hamiltonian is

Ĥ(x) = Ĥ0 + V (x), (1)
and the total energy can be expanded with respect to x as

E(x) = Ψ(x)| Ĥ(x)|Ψ(x) Ψ(x)|Ψ(x) = E (0) + E (1) x + 1 2 E (2) x 2 + • • • , (2)
where Ψ(x) is the wave function of the state considered, which is usually the ground state. In this expansion, E (0) = E 0 is the unperturbed energy associated to the Hamiltonian Ĥ0 , and the energy derivatives

E (1) = dE dx x=0 ; E (2) = d 2 E dx 2 x=0 ; etc... (3)
are called time-independent molecular properties. They are characteristic of the molecule and its quantum state, and contain important information about the response of the system to the perturbation.

Examples 1.2.1 Geometrical derivatives

One the most important example of a static perturbation is a deformation of the molecular geometry (in the Born-Oppenheimer approximation). In this case, x corresponds to variations of nuclei positions ∆R i , and the total energy is expanded as

E({R i }) = E (0) + i E (1)
i ∆R i + 1 2 i,j E (2)
i,j

∆R i ∆R j + • • • , (4)
where

E (1) i = ∂E ∂R i (5)
is called the molecular gradient, and

E (2) i,j = ∂ 2 E ∂R i ∂R j (6)
is called the molecular Hessian. These quantities are used in geometry optimization for locating and characterizing critical points on the molecular potential energy surface (energy minimum, saddle point at a transition state). There are also used for calculating spectroscopic constants such as harmonic vibrational frequencies. Note that anharmonicity vibrational corrections require the calculation of third-order derivatives E

i,j,k = ∂ 3 E/∂R i ∂R j ∂R k .

Electric properties

In this case, the perturbation applied is an external electrostatic potential v(r). Let us consider the corresponding interaction energy

E int = n(r)v(r)dr, (7)
where n(r) is the charge density of the system. Expanding the potential v(r) around an origin point r = 0, we have

v(r) = v(0) + i ∂v ∂r i r=0 r i + 1 2 i,j ∂ 2 v ∂r i ∂r j r=0 r i r j + • • • = v(0) - i E i r i - 1 2 i,j F i,j r i r j + • • • , (8)
where v(0) is a constant potential, E i are the components of the corresponding electric field, and F i,j are the components of the electric field gradient. The interaction energy can thus be expanded as

E int = q v(0) - i µ i E i - 1 2 i,j Q i,j F i,j + • • • , (9)
where q = n(r)dr is the electric charge of the system, µ i = n(r)r i dr are the components of the electric dipole moment, and Q i,j = n(r)r i r j dr are the components of electric quadrupole moment. Hence, these properties can be defined as energy derivatives q = ∂E ∂v(0)

;

µ i = - ∂E ∂E i ; Q i,j = -2 ∂E ∂F i,j . (10)
We can further decompose the dipole moment into permanent and field-induced contributions

µ i = µ 0,i + j α i,j E j + • • • , (11)
which leads to the definition of the permanent electric dipole moment

µ 0,i = - ∂E ∂E i E=0 , (12)
and the electric dipole polarizability

α i,j = - ∂ 2 E ∂E i ∂E j E=0 . ( 13
)

Magnetic properties

For magnetic properties, we consider the double expansion of the total energy with respect to nuclear magnetic dipole moments m a (a refers to an atomic nucleus) and an externally applied magnetic field B

E({m a } , B) = E (0) + a E (10) a m a + E (01) B + 1 2 a,b E (20) a,b m a m b + a E (11) a m a B + 1 2 E (02) B 2 + • • • , (14)
the dependence of the total energy on the nuclear magnetic dipole moments requiring the inclusion of the hyperfine-interaction terms in the Hamiltonian.

In nuclear magnetic resonance (NMR), we measure the coupling between m a and B, which is quantified by the second-order derivative

E (11) a = ∂ 2 E ∂m a ∂B . (15)
In vacuum, the interaction energy is simply

E vac = -m a B, giving E (11)
a,vac = -1. In a molecule, due to the presence of electrons, it is modified by a few ppm, and we write E (11) a,mol = -1 + σ a where σ a is the NMR shielding constant of the nucleus a. We also measure the nuclear spin-spin coupling constants between two nuclei

E (20) a,b = ∂ 2 E ∂m a ∂m b = J a,b . (16)
In solution, only the indirect contribution to this coupling (mediated by the electrons) survives and provides useful geometrical information about the molecule.

Similarly to the corresponding electric properties, we can also define properties such as the permanent magnetic dipole moment M 0 = -∂E/∂B| B=0 and the dipole magnetizability

ξ = -∂ 2 E/∂B 2 | B=0 .
There are also other important magnetic properties used in electron paramagnetic resonance (EPR) such as EPR hyperfine coupling constants and g values.

Calculation of energy derivatives 1.3.1 Numerical versus analytical differentiation

There are two ways of calculating derivatives, by numerical differentiation and by analytical differentiation.

In numerical differentiation, the derivatives are calculated by finite differences or by polynomial fitting.

• It is simple to implement.

• The numerical precision is limited.

• The computational efficiency is low.

• It is not general. Usually, only real-valued and static perturbations can be done.

In analytical differentiation, the derivatives are calculated explicitly from the analytical expressions.

• It is difficult to implement.

• The precision is higher.

• The computational efficiency is higher.

• It is more general. In particular, frequency-dependent perturbations are possible.

Overall, analytical differentiation is preferable.

Sum-over-state expression for exact wave functions

For "exact" wave functions (full configuration interaction in a basis), we can use straightforward perturbation theory to find the expression of the energy derivatives. Expanding the perturbation operator V (x) in powers of x, we have

Ĥ(x) = Ĥ0 + V (1) x + 1 2 V (2) x 2 + • • • . ( 17
)
The first-order energy derivative is

E (1) = Ψ 0 | V (1) |Ψ 0 , (18)
where Ψ 0 is the exact wave function of the unperturbed Hamiltonian Ĥ0 for the state considered, normalized such that Ψ 0 |Ψ 0 = 1. Similarly, the second-order energy derivative is found from the first-order wave function correction Ψ (1) , using intermediate normalization Ψ 0 |Ψ(x) = 1,

E (2) = Ψ 0 | V (2) |Ψ 0 + 2 Ψ 0 | V (1) |Ψ (1) = Ψ 0 | V (2) |Ψ 0 -2 n =0 Ψ 0 | V (1) |Ψ n Ψ n | V (1) |Ψ 0 E n -E 0 , (19)
where Ψ n and E n form a set of complete exact eigenstates and associated eigenvalues of Ĥ0 . This last expression is very impractical since we need to know all the exact states of Ĥ0 , which is a complicated many-body Hamiltonian. Moreover, for approximate methods such as Hartree-Fock, it is not clear what we should use for Ψ n and E n .

General expressions for approximate methods

In an (approximate) electronic-structure method, the total energy E(x) is obtained by optimizing parameters p = (p 1 , p 2 , ....) in an energy function E(x, p) for each fixed value of x. The final total energy E(x) is obtained for the optimal value of the parameters p o (x) which are functions of x

E(x) = E(x, p o (x)). ( 20
)
Note that the optimization is not necessarily variational, i.e. the optimization criterion is not necessarily to minimize E(x, p) with respect to the parameters p.

Examples:

• Hartree-Fock (HF), Kohn-Sham (KS): p are the orbital parameters. The parameters are variational.

• Multiconfiguration self-consistent field (MCSCF): p are the coefficients of the configurations (or determinants) and the orbital parameters. All parameters are variational.

• Configuration interaction (CI): p are the coefficients of the configurations (or determinants) and the orbital parameters. The coefficients of the configurations are variational but not the orbital parameters.

• Coupled cluster (CC): p are the cluster amplitudes and the orbital parameters. None parameters are variational.

First-order energy derivative

The first-order derivative of E(x) with respect to x is made of two terms

dE(x) dx = ∂E(x, p o ) ∂x + i ∂E(x, p) ∂p i p=p o ∂p o i ∂x , (21)
where the first term constitutes the explicit dependence on x and the second term constitutes the implicit dependence on x. The derivative ∂p o i /∂x appearing in the second term is called the wave-function linear-response vector and contains information about how the electronic structure changes when the system is perturbed. It is not straightforward to calculate it since we do not know explicitly dependence of p o i on x. Where all the parameters are variational, there is an important simplification. Indeed, in this case, the stationary condition (zero electronic gradient)

∂E(x, p) ∂p i p=p o = 0, (22)
implies that the second term in Eq. ( 21) vanishes, so the first-order energy derivative reduces to

dE(x) dx = ∂E(x, p o ) ∂x , (23)
and we do not need the wave-function linear-response vector. For example, for the calculation of the HF molecular gradient, we need to consider only the explicit dependence on the nuclear positions in the Hamiltonian and in the wave function, but not the implicit dependence of the orbital coefficients on the nuclear positions since these parameters are variational.

Second-order energy derivative

We will only consider the case where all the parameters are variational. In this case, the secondorder derivative of E(x) is obtained by taking the derivative of Eq. (23)

d 2 E(x) dx 2 = ∂ 2 E(x, p o ) ∂x 2 + i ∂ 2 E(x, p) ∂x∂p i p=p o ∂p o i ∂x , (24)
where the perturbed electronic gradient ∂ 2 E(x, p)/∂x∂p i | p=p o is not zero. We thus now need the wave-function linear-response vector. Note that we need only the first-order derivative of the wave function to calculate the secondorder derivative of the energy. More generally, we have the so-called 2n + 1 rule: For variational parameters, the derivatives of the wave function to order n determine the energy derivatives to order 2n + 1.

Linear response equations

In order to obtain the wave-function linear-response vector ∂p o i /∂x in the case where all the parameters are variational, we start by noting that the stationary condition is true for all x ∀x, ∂E(x, p)

∂p i p=p o = 0, (25)
which means that we can take the first-order derivative of Eq. ( 25) with respect to x, giving

∂ 2 E(x, p) ∂x∂p i p=p o + j ∂ 2 E(x, p) ∂p i ∂p j p=p o ∂p o j ∂x = 0. ( 26
)
We thus arrive at the linear response equations

j ∂ 2 E(x, p) ∂p i ∂p j p=p o ∂p o j ∂x = - ∂ 2 E(x, p) ∂x∂p i p=p o , (27)
which is a linear system of equations whose solution gives the wave-function linear-response vector ∂p o i /∂x. Remember that, in practice, we are interested in energy derivatives evaluated at x = 0. Therefore, these equations must in fact be solved using the unperturbed electronic Hessian ∂ 2 E(x = 0, p)/∂p i ∂p j (independent of the perturbation). Since its dimensions are usually large, the response equations are usually solved iteratively without explicitly constructing and storing the Hessian matrix.

Remark: Lagrangian formalism for non-variational parameters

For methods with non-variational parameters (CI, CC, MP2, ...), we can simplify the calculation by using the technique of Lagrange's multipliers. Let us consider the case of CI with parameters p = (c, κ) where c are the coefficients of the configurations which are variational, and κ are the orbital rotation parameters which are non variational. Indeed, the orbital rotation parameters makes the HF energy stationary, ∂E HF (x, κ)/∂κ i = 0, but not the CI energy, ∂E CI (x, c, κ)/∂κ i = 0. In this case, we define the Lagrange function L CI (x, c, κ, κ) depending on Lagrange's multipliers κ

L CI (x, c, κ, κ) = E CI (x, c, κ) + i κi ∂E HF (x, κ) ∂κ i . (28)
For all values of x, L CI (x, c, κ, κ) takes the same value than E CI (x, c, κ), but L CI (x, c, κ, κ) has the advantage of being stationary with respect to all parameters c, κ and κ. This simplifies the calculation of the derivatives with respect to x, at the price of determining the Lagrange's multipliers κ.

1.4 Example of linear-response equations for Hartree-Fock: Coupled-perturbed Hartree-Fock

Exponential parametrization

For optimizing the orbitals in the HF determinant wave function, it is convenient to use an exponential parametrization

|Φ(κ) = e κ|Φ 0 , (29)
where e κ is an unitary operator performing rotations between occupied and virtual spin orbitals in a reference determinant wave function Φ 0 . This rotation operator is constructed from an anti-Hermitian single-excitation operator, κ = -κ † , which can be written in the second-quantization formalism as

κ = occ a vir r κ ar â † r âa -κ * ar â † a âr , (30)
where â † k and âk are the creation and annihilation operators of the spin orbital k, respectively, and the indices a and r refer to occupied and virtual spin orbitals in the reference determinant, respectively. The parameters κ = {κ ar } are called the orbital rotation parameters. In comparison to the orbital coefficients on the atomic basis functions, the orbital rotation parameters have the advantage of providing a non-redundant parametrization of the wave function, so that one can vary them independently without having to impose any constraints.

At each step of the optimization, the orbitals in the reference determinant Φ 0 are updated so that one always considers variations of the orbital rotation parameters around κ = 0. The expansion of the HF wave function with respect to κ thus writes

|Φ(κ) = 1 + κ + 1 2 κ2 + • • • |Φ 0 = |Φ 0 +
where |Φ r a = â † r âa |Φ 0 is a singly excited determinant and |Φ rs ab = â † r âa â † s âb |Φ 0 is a doubly excited determinant. Note that |Φ rs ab = 0 if a = b or r = s, since a given spin orbital cannot be annihilated or created more than once.

Electronic gradient and Hessian

We now want to calculate the HF electronic gradient and Hessian, i.e. the first-and second-order derivatives of the HF total energy

E HF (κ) = Φ(κ)| Ĥ0 |Φ(κ) Φ(κ)|Φ(κ) (32)
with respect to the parameters κ. Here, we will only consider the case of real-valued orbitals and real-valued orbital parameters κ * ar = κ ar . In this case, from Eq. ( 31), the first-and second-order derivatives of the HF wave function |Φ(κ) are

∂|Φ(κ) ∂κ ar κ=0 = |Φ r a , (33)
and

∂ 2 |Φ(κ) ∂κ ar ∂κ bs κ=0 = | Φrs ab = |Φ rs ab if a = b or r = s -|Φ 0 if a = b and r = s. ( 34
)
Exercise: Assuming the reference wave function Φ 0 is normalized, Φ 0 |Φ 0 = 1, show that the HF electronic gradient writes

∂E HF ∂κ ar κ=0 = 2 Φ r a | Ĥ0 |Φ 0 , (35)
and is thus zero for converged HF spin orbitals (Brillouin's theorem). Also show that the HF electronic Hessian writes

∂
A ar,bs = (ε r -ε a ) δ ab δ rs + rb||as , (37)
and

B ar,bs = rs||ab , (38)
where ε k is the energy of the HF spin orbital k, and rb||as = rb|as -rb|sa are the antisymmetrized two-electron integrals over the HF spin orbitals.

Example: Calculation of the electric dipole polarizability

The so-called coupled-perturbed Hartree-Fock (CPHF) electric dipole polarizability is, according to Eq. ( 24),

α CPHF i,j = - ∂ 2 E HF ∂E i ∂E j E=0 = - ∂ 2 E HF ∂E i ∂E j + occ a vir r ∂ 2 E HF ∂E i ∂κ ar κ=0 ∂κ ar ∂E j , (39)
where in this expression E HF is the HF energy in the presence of a perturbing uniform electric field E

E HF (E, κ) = Φ(κ)| Ĥ0 -μ • E|Φ(κ) Φ(κ)|Φ(κ) , ( 40
)
where μ is the dipole moment operator.

Clearly, since the explicit dependence of HF energy on the electric field is linear, we have ∂ 2 E HF /∂E i ∂E j = 0. Similarly to the calculation leading to Eq. ( 35), the HF perturbed electronic gradient is found to be

∂ 2 E HF ∂E i ∂κ ar κ=0 = -2 Φ r a |μ i |Φ 0 = -2 r|μ i |a , (41)
where r|μ i |a are the transition dipole moment one-electron integrals. Therefore, Eq. (39) gives

α CPHF i,j = occ a vir r 2 r|μ i |a ∂κ ar ∂E j , (42)
and it only remains to calculate the wave-function linear-response vector ∂κ ar /∂E j . This requires the resolution of the HF linear-response equations, also known as the CPHF equations, which write, using Eqs. ( 27), ( 36) and (41),

Using vector/matrix notations, the CPHF electric dipole polarizability finally writes

α CPHF i,j = 2µ T i (A + B) -1 µ j , (44)
where µ i is the vector of components µ i,ar = r|μ i |a and A + B is the matrix of elements A ar,bs + B ar,bs .

Remark: The non-interacting case: Uncoupled Hartree-Fock

In the simplified case of no electron-electron interaction, also called uncoupled Hartree-Fock (UCHF) in this context, we have A ar,bs = (ε r -ε a ) δ ab δ rs and B ar,bs = 0. Therefore, Eq. ( 44) becomes

α UCHF i,j = 2 occ a vir r r|μ i |a a|μ j |r ε r -ε a . ( 45
)
This simple expression is convenient for discussing the physical contents of the polarizability. Indeed, it is clearly seen that the UCHF static dipole polarizability is proportional to the transition dipole-moment integrals r|μ i |a and inversely proportional to the (non-interacting) excitation energies ε r -ε a . Therefore, we expect the polarizability to be large for systems with low-lying excited states that can be reached from the ground-state by dipole-allowed transitions. For example, alkaline-earth atoms (Be, Mg, Ca, ...) with their low-lying 2p orbitals have much larger polarizabilities than the rare-gas atoms (He, Ne, Ar, Kr, ...) which have large gaps.

Examples of spherically averaged static dipole polarizabilities ᾱ = (α x,x + α y,y + α z,z )/3 for rare-gas and alkaline-earth-metal atoms calculated at the UCHF and CPHF levels are given in Table 1. The polarizabilities are indeed one order of magnitude larger for alkaline-earth-metal atoms than for rare-gas atoms. The CPHF polarizabilities are in relatively good agreement with the estimated exact values for the rare-gas atoms. The UCHF or CPHF polarizabilities are less accurate for alkaline-earth-metal atoms. This is due to the neglect of electron correlation effects, which are particularly large in these systems because of their small HF-orbital gaps. 2 Time-dependent molecular properties

Time-dependent Schrödinger equation and quasi-energy

We now consider a time-dependent Hamiltonian with a periodic time-dependent perturbation operator of the form

Ĥ(t) = Ĥ0 + x 1 V1 e -iωt + x 2 V2 e +iωt , ( 46
)
where Ĥ0 is still the time-independent unperturbed Hamiltonian, x 1 and x 2 are (real-valued) variables representing the strengths of each perturbation term, and ω is the pulsation (or frequency). For the perturbation operator to be Hermitian, we must impose the condition V1 = V † 2 and note that we will ultimately be interested in the case x 1 = x 2 . Think of this perturbation as representing the interaction of a molecule with a monochromatic electromagnetic field.

The time-dependent wave function Ψ(t) of this system satisfies the time-dependent Schrödinger equation

Ĥ(t) -i ∂ ∂t | Ψ(t) = 0. ( 47
)
This equation can be reformulated in a more convenient form. We first extract a phase factor F(t) out of the wave function

| Ψ(t) = e -iF (t) |Ψ(t) , (48)
where Ψ(t) is a new time-dependent wave function which satisfies the equation

Ĥ(t) -i ∂ ∂t |Ψ(t) = Ḟ (t)|Ψ(t) , (49)
where Ḟ (t) = dF/dt. Equation ( 49) is in the form of an eigenvalue equation, similar to the time-independent Schrödinger equation, and Ḟ (t) plays a role similar to the energy. As the Hamiltonian is time periodic, Ĥ(t + T ) = Ĥ(t) where T = 2π/ω is the period, we define the quasi-energy as the time average of Ḟ(t) over a period

Q = 1 T T 0 dt Ḟ(t) = 1 T T 0 dt Ψ(t)| Ĥ(t) -i ∂ ∂t |Ψ(t) Ψ(t)|Ψ(t) . ( 50
)
The advantage of introducing the quasi-energy is that there is variational principle similar to the time-independent one, namely the exact wave function Ψ(t) solution of Eq. (49) makes Q stationary with respect to variations in Ψ(t). In the time-independent case, the quasi-energy Q reduces to the energy E and the wave function Ψ(t) reduces to the time-independent wave function Ψ.

Definition of time-dependent molecular properties

The quasi-energy can be expanded with respect to the perturbation parameters x 1 and x 2

Q = Q (0) + Q (10) x 1 + Q (01) x 2 + 1 2 Q (20) x 2 1 + Q (11) x 1 x 2 + 1 2 Q (02) x 2 2 + • • • , (51)
where

Q (0) = Ψ 0 | Ĥ0 |Ψ 0 = E (0) = E 0
is the unperturbed energy corresponding to the unperturbed normalized state considered Ψ 0 , and the derivatives

Q (10) = ∂Q ∂x 1 x=0 ; Q (01) = ∂Q ∂x 2 x=0 Q (20) = ∂ 2 Q ∂x 2 1 x=0 ; Q (11) = 2 ∂ 2 Q ∂x 1 ∂x 2 x=0 ; Q (02) = ∂ 2 Q ∂x 2 2 x=0
; etc..,

define time-dependent molecular properties. They depend on the frequency ω of the perturbation, and encompass the static molecular properties as a special case for ω = 0.

Sum-over-state expressions for exact wave functions

Using standard time-dependent perturbation theory on Eq. ( 49), we obtain the first-order derivative Q (10)

Q (10) = 1 T T 0 dt Ψ 0 | V1 e -iωt |Ψ 0 = Ψ 0 | V1 |Ψ 0 δ(ω), (53)
where we used (1/T ) T 0 dt e -iωt = δ(ω). Here, δ(ω) must be understood as a Kronecker delta for finite T and as a Dirac delta function for T → ∞. It is seen that only the static component (ω = 0) is not zero for first-order molecular properties. This can be intuitively understood as follows: We consider molecular properties for a fixed state |Ψ 0 , and thus "absorption" or "emission" of a photon of frequency ω = 0 is not allowed by energy conservation.

Similarly, we find the second-order derivative Q (20) by perturbation theory on Eq. ( 49), using intermediate normalization Ψ 0 |Ψ(t) = 1,

Q (20) = 1 T T 0 dt 2 Ψ 0 | V1 e -iωt |Ψ (10) (t) , (54)
where the first-order wave function correction has the form

|Ψ (10) (t) = n =0 c (10) n (t)e -iωnt |Ψ n , (55)
where ω n = E n -E 0 are the excitation energies of the unperturbed system, and the first-order coefficients c

(10) n (t) are found using Eq. ( 49)

c (10) n (t) = -i t -∞ dt ′ Ψ n | V1 e -iωt ′ e γt ′ |Ψ 0 e iωnt ′ , (56)
where the factor e γt ′ with γ → 0 + adiabatically switches on the perturbation from t ′ → -∞, imposing the initial condition c

n (t) = Ψ n | V1 e -iωt |Ψ 0 e iωnt ω -ω n + iγ , (57)
where the factor e γt has been dropped. Combining Eqs. (54), ( 55) and (57), we find

Q (20) = 1 T T 0 dt 2 n =0 Ψ 0 | V1 e -iωt |Ψ n Ψ n | V1 e -iωt |Ψ 0 ω -ω n + iγ = 2 n =0 Ψ 0 | V1 |Ψ n Ψ n | V1 |Ψ 0 ω -ω n + iγ δ(2ω), (58)
where again only the static component is not zero. The case of the mixed second-order derivative Q (11) is more interesting. The calculation is similar to the previous one and gives

Q (11) = 1 T T 0 dt Ψ 0 | V2 e +iωt |Ψ (10) (t) + Ψ 0 | V1 e -iωt |Ψ (01) (t) = 1 T T 0 dt n =0 Ψ 0 | V2 e +iωt |Ψ n Ψ n | V1 e -iωt |Ψ 0 ω -ω n + iγ + n =0 Ψ 0 | V1 e -iωt |Ψ n Ψ n | V2 e +iωt |Ψ 0 -ω -ω n + iγ = n =0 Ψ 0 | V2 |Ψ n Ψ n | V1 |Ψ 0 ω -ω n + iγ - n =0 Ψ 0 | V1 |Ψ n Ψ n | V2 |Ψ 0 ω + ω n -iγ , ( 59
)
where we see that Q (11) does not vanish for ω = 0. This is truly a dynamic molecular property. Physically, it corresponds to the following phenomenon: The "absorption" of a (virtual) photon of frequency ω (due to the perturbation V1 e -iωt ) followed by the "emission" of a (virtual) photon of the same frequency (due to the perturbation V2 e +iωt ), which is allowed since there is no net change in energy, probing in the process the excited states of the system. The quantity Q (11) is called the linear-response function of operators V2 and V1 in the spectral representation and it is often denoted as Q (11) = V2 ; V1 ω . It can also be seen as giving the change in the expectation value of the operator V2 due to the perturbation V1 e -iωt . The linear-response function has poles in ω at the excitation energies of the unperturbed system ω n = E n -E 0 (and at the deexcitation energies -ω n for the second term), and hence gives access to all the excitation energies of the system. The infinitesimal imaginary shifts ±iγ are kept because they are important for manipulation of the linear-response function using complex analysis. Similarly, higher-order derivatives of the quasi-energy Q leads to quadratic-response functions, cubic-response functions, etc...

Example of the dynamic electric dipole polarizability

The simplest example of a second-order dynamic molecular property is the dynamic electric dipole polarizability. We consider the following perturbed Hamiltonian

Ĥ(t) = Ĥ0 -E + • μ e -iωt -E -• μ e +iωt , ( 60
)
where E + and E -are the amplitudes of the oscillating electric field and μ is the dipole-moment operator. The amplitudes E + and E -are distinguished for the theoretical derivation but we will eventually be interested in the physical case E + = E -= E. The dynamic electric dipole polarizability is then defined as the mixed second-order derivative

α i,j (ω) = - ∂ 2 Q ∂E - i ∂E + j E=0 , (61)
where E = 0 refers to E + = E -= 0. Using Eq. ( 59) with V2 = -μ i and V1 = -μ j , we find

α i,j (ω) = - n =0 Ψ 0 |μ i |Ψ n Ψ n |μ j |Ψ 0 ω -ω n + iγ + n =0 Ψ 0 |μ j |Ψ n Ψ n |μ i |Ψ 0 ω + ω n -iγ . ( 62
)
We are often interested only in the spherical average (or isotropic component) of the polarizability tensor ᾱ(ω) = (1/3) i α i,i (ω), which can be written as (disregarding the infinitesimal shifts ±iγ)

ᾱ(ω) = n =0 f n ω 2 n -ω 2 , ( 63
)
where f n are the oscillator strengths

f n = 2 3 ω n i | Ψ 0 |μ i |Ψ n | 2 . ( 64
)
The oscillator strength f n gives the intensity in the absorption or emission spectrum of the electronic transition from state |Ψ 0 to state |Ψ n corresponding to the excitation energy ω n . The knowledge of the dynamic polarizability ᾱ(ω) thus gives access to the full absorption or emission spectrum of the system.

General expressions for approximate methods

Just like in the static case, the previous expressions of the dynamic molecular properties are not applicable for approximate methods such HF, MCSCF, etc... Due to the analogy between the quasi-energy Q in the dynamic case and the energy E in the static case, essentially the same expressions derived in Section 1.3.3 for approximate methods apply here as well.

The quasi-energy Q is obtained by plugging optimal parameters p o (x 1 , x 2 ) in a quasi-energy function Q(x 1 , x 2 ; p)

Q(x 1 , x 2 ) = Q(x 1 , x 2 ; p o (x 1 , x 2 )). (65)
In contrast to the static case, the parameters now depend on time and take complex values. They are written as

p = c e -iωt + d * e +iωt , (66)
where c and d * are complex numbers. The complex conjugate d * is introduced instead of d for symmetry reason. We will consider only the case where all the parameters are variational, i.e. fulfilling the stationary conditions

∂Q ∂c * i p=p o = 0 and ∂Q ∂d * i p=p o = 0, (67)
and similarly for the complex-conjugate equations. Note that the case of non-variational parameters can be treated using the Lagrangian formalism. Thus, only the explicit dependence term contributes to the first-order derivative of the quasi-energy

∂Q ∂x 1 = ∂Q ∂x 1 , (68)
and the second-order derivative of the quasi-energy is

∂ 2 Q ∂x 1 ∂x 2 = ∂ 2 Q ∂x 1 ∂x 2 + i ∂ 2 Q ∂x 1 ∂c i ∂c i ∂x 2 + ∂ 2 Q ∂x 1 ∂d i ∂d i ∂x 2 + ∂ 2 Q ∂x 1 ∂c * i ∂c * i ∂x 2 + ∂ 2 Q ∂x 1 ∂d * i ∂d * i ∂x 2 , ( 69
)
where it is understood that the parameters are evaluated at their optimal values p = p o . The wave-function linear-response vectors ∂c i /∂x 2 and ∂d i /∂x 2 are determined from the timedependent linear-response equations, written in a 2 × 2-matrix form,

j   ∂ 2 Q ∂c * i ∂c j ∂ 2 Q ∂c * i ∂d j ∂ 2 Q ∂d * i ∂c j ∂ 2 Q ∂d * i ∂d j     ∂c j ∂x 2 ∂d j ∂x 2   = -   ∂ 2 Q ∂x 2 ∂c * i ∂ 2 Q ∂x 2 ∂d * i   , (70)
where it has been assumed that the terms

∂ 2 Q/∂c * i ∂c * j , ∂ 2 Q/∂d * i ∂d * j , and ∂ 2 Q/∂c * i ∂d * j vanish
. This is indeed normally the case, as it will be illustrated in the next section for time-dependent Hartree-Fock. The linear-response vectors ∂c * i /∂x 2 and ∂d * i /∂x 2 satisfy a similar equation.

2.5

Example: Time-dependent Hartree-Fock

Time-dependent electronic Hessian

According to Eq. (66), the orbital rotation parameters are written as

κ = c e -iωt + d * e +iωt , (71)
and we want to calculate the second-order derivatives of the HF quasi-energy

Q HF (κ) = 1 T T 0 dt Φ(κ)| Ĥ0 -i ∂ ∂t |Φ(κ) Φ(κ)|Φ(κ) , (72)
with respect to the parameters c, c * , d, and d * . We will work out the general case of complexvalued HF orbitals. By plugging the time-dependent parameters of Eq. ( 71) in the expansion of the HF wave function |Φ(κ) in Eq. (31), we find that the only non-zero first-and second-order derivatives of |Φ(κ) are ∂|Φ(κ)

∂c ar κ=0 = e -iωt |Φ r a , (73)
∂|Φ(κ) ∂d * ar κ=0 = e +iωt |Φ r a , (74)
and

∂ 2 |Φ(κ) ∂c ar ∂c bs κ=0 = e -2iωt |Φ rs ab , (75)
∂ 2 |Φ(κ) ∂d * ar ∂d * bs κ=0 = e +2iωt |Φ rs ab , (76)
∂ 2 |Φ(κ) ∂c ar ∂d * bs κ=0 = ∂ 2 |Φ(κ) ∂d * ar ∂c bs κ=0 = |Φ rs ab , (77)
∂ 2 |Φ(κ) ∂c * ar ∂c ar κ=0 = ∂ 2 |Φ(κ) ∂d * ar ∂d ar κ=0 = -|Φ 0 , (78)
∂ 2 |Φ(κ) ∂c ar ∂d ar κ=0 = -e -2iωt |Φ 0 , (79)
∂ 2 |Φ(κ) ∂c * ar ∂d * ar κ=0 = -e +2iωt |Φ 0 . (80)
Exercise: By using Eqs. ( 73)-(80) and performing the time average in Eq. ( 72), show that the second-order derivatives of the HF quasi-energy are, for ω = 0,

∂ 2 Q HF ∂c * ar ∂c bs κ=0 = A ar,bs -ωS ar,bs , (81)
∂ 2 Q HF ∂d * ar ∂d bs κ=0 = A * ar,bs + ωS * ar,bs , (82)
∂ 2 Q HF ∂c * ar ∂d bs κ=0 = B ar,bs , (83)
∂ 2 Q HF ∂d * ar ∂c bs κ=0 = B * ar,bs , (84)

Time-dependent linear-response equation

In matrix/vector notations, the time-dependent linear-response equation of Eq. (70) takes the following form in HF

A B B * A * -ω 1 0 0 -1 X Y = - V W , (85)
where A and B are the matrices of elements A ar,bs and B ar,bs , respectively, 1 is the identity matrix, (X, Y) is the linear-response vector and (V, W) is the perturbed gradient vector. It is interesting to consider the linear-response equation without the right-hand-side perturbation term, which gives the following generalized eigenvalue equation

A B B * A * X n Y n = ω n 1 0 0 -1 X n Y n . ( 86
)
Note that the matrix A is Hermitian and the matrix B is symmetric, and thus the 2 × 2-block matrix on the left-hand-side of Eq. ( 86) is Hermitian. However, the 2 × 2-block matrix on the right-hand-side of Eq. ( 86) is not positive definite. Eq. ( 86) is not a Hermitian eigenvalue equation but a pseudo-Hermitian eigenvalue equation. One can show that if the 2 × 2-block matrix on the left-hand-side of Eq. ( 86) is positive definite then all the eigenvalues ω n are real numbers. Since this matrix corresponds to the Hessian of the HF energy with respect to the orbital rotation parameters, it is positive definite if the HF solution is a stable minimum with respect to the variation of all orbital parameters. This is known as the Hartree-Fock stability conditions.

The solutions of Eq. ( 86) come in pairs: If (X n , Y n ) is an eigenvector with eigenvalue ω n , then it is easy to check that (Y * n , X * n ) is also an eigenvector with opposite eigenvalue -ω n . When calculating the dynamic dipole polarizability, we will see that the positive eigenvalues, ω n > 0, correspond to excitation energies in the time-dependent Hartree-Fock (TDHF) approximation, while the negative eigenvalues, ω n < 0, correspond to de-excitation energies.

Dynamic dipole polarizability

As a second-order derivative, the TDHF dynamic dipole polarizability is calculated according to Eq. ( 69)

α TDHF i,j (ω) = - ∂ 2 Q HF ∂E - i ∂E + j = - ∂ 2 Q HF ∂E - i ∂E + j + occ a vir r ∂ 2 Q HF ∂E - i ∂c ar ∂c ar ∂E + j + ∂ 2 Q HF ∂E - i ∂d ar ∂d ar ∂E + j + occ a vir r ∂ 2 Q HF ∂E - i ∂c * ar ∂c * ar ∂E + j + ∂ 2 Q HF ∂E - i ∂d * ar ∂d * ar ∂E + j , (87)
where all the terms are evaluated at E = 0 and κ = 0. In this expression Q HF is the HF quasienergy of the perturbed system

Q HF (E, κ) = 1 T T 0 dt Φ(κ)| Ĥ(t) -i ∂ ∂t |Φ(κ) Φ(κ)|Φ(κ) , (88)
with the Hamiltonian Ĥ(t) in Eq. (60) which includes the oscillating electric field. As for the static case, the second-order derivative ∂ 2 Q HF /∂E - i ∂E + j vanishes, and the first two perturbed gradients in Eq. (87) give the transition dipole-moment integrals

∂ 2 Q HF ∂E - i ∂c ar = -Φ 0 |μ i |Φ r a = -a|μ i |r , (89)
∂ 2 Q HF ∂E - i ∂d ar = -Φ r a |μ i |Φ 0 = -r|μ i |a , (90)
the two other terms being zero, ∂ 2 Q HF /∂E - i ∂c * ar = ∂ 2 Q HF /∂E - i ∂d * ar = 0. This leads to the following expression for the polarizability

α TDHF i,j (ω) = µ i µ * i † A B B * A * -ω 1 0 0 -1 -1 µ j µ * j , (91)
where µ i is the vector of elements µ i,ar = r|μ i |a .

Exercise: Check that the inverse of the time-dependent Hessian can be written as the following spectral representation:

A B B * A * -ω 1 0 0 -1 -1 = n 1 ω n -ω X n Y n X n Y n † + 1 ω n + ω Y * n X * n Y * n X * n † , (92)
where the sum is over eigenvectors with positive eigenvalues ω n > 0 only. The following completeness relation will be assumed:

1 0 0 1 = n X n -Y n X n Y n † - Y * n -X * n Y * n X * n † . ( 93
)
Plugging the spectral representation of Eq. (92) into Eq. ( 91), we finally arrive at, for realvalued HF orbitals,

α TDHF i,j (ω) = - n µ T i (X n + Y n ) (X n + Y n ) T µ j ω -ω n + µ T i (X n + Y n ) (X n + Y n ) T µ j ω + ω n . ( 94
)
The comparison with the exact expression of the dynamic polarizability in Eq. (62) then confirms that the eigenvalues ω n in Eq. (86) correspond to excitation energies. The corresponding TDHF oscillator strengths are

f n = 2 3 ω n i (X n + Y n ) T µ i 2 .
(95)

Photoabsorption spectrum

TDHF can be used to calculate the discrete part of the photoabsorption spectrum of an electronic system. Indeed, the resolution of the linear-response equation, Eq. ( 86), directly leads to the excitation energies ω n from the ground state to state n, and to the associated oscillator strengths f n (related to the probability of the transitions and therefore to the intensity of the spectrum). Note, however, that linear-response theory does not provide the wave functions of the excited states. This leads to the question of how to assign the calculated excitation energies to a physical state.

It is useful to consider the so-called Tamm-Dancoff approximation (TDA) to linear-response theory which corresponds to neglecting the matrix B in Eq. (86). After setting B = 0, Eq. (86) simplifies to a standard Hermitian eigenvalue equation

AX n = ω n X n , (96)
where we recall that A is the matrix of the Hamiltonian (shifted by the HF energy) in the basis of the single-excited determinants: A ar,bs = Φ r a | Ĥ0 -E HF |Φ s b . Eq. ( 96) is simply a configuration-interaction calculation with only single excitations (CIS), whose resolution gives the excited-state wave functions expanded in the basis of single-excited determinants

|Ψ n = occ a vir r X n,ar |Φ r a , (97)
and X n,ar is the coefficient for the single orbital excitation a → r. Thus, in this case, the coefficients X n exactly provide the physical contents of the excited state n as contributions of single orbital excitations. Usually, one orbital excitation dominates over all the other orbital excitations, and this is used to assign the state to this orbital excitation. By analogy, in the general linear-response theory case without the TDA, X n,ar + Y n,ar is often taken as an approximation to the coefficient of the wave function of the excited state n on the single-excited determinant Φ r a . Table 2 shows an example of a TDHF calculation of excitation energies and oscillator strengths, without and with the TDA, on the CO molecule. The results are compared with those obtained with time-dependent density-functional theory (TDDFT) with the local-density approximation (LDA), which shares essentially the same formalism, and with a more accurate equation-of-motion coupled-cluster singles doubles (EOM-CCSD) calculation. Experimental values are also given. Due to the difficulty of extracting vertical excitation energies and oscillator strengths from the experimental spectrum, the EOM-CCSD values are taken as reference.

The TDHF excitation energies to the low-lying valence spin-triplet states are very much underestimated, while the TDHF excitation energies to the valence spin-singlet states tend to be more accurate. The TDHF excitation energies to the (diffuse and high-lying) Rydberg states tend to be overestimated. The TDHF underestimation of the excitation energies to the valence spin-triplet states is connected to the so-called near-triplet Hartree-Fock instability. This refers to the fact that, for this system, the spin-restricted Hartree-Fock (RHF) ground-state wave function is close to being unstable to spin-triplet variations of the orbitals breaking the spin symmetry. The TDA largely reduces the TDHF underestimation of the excitation energies to the valence spin-triplet states, while not changing much the excitation energies to Rydberg states.

TDLDA does not suffer from near-triplet Hartree-Fock instability and overall gives more accurate valence excitation energies than TDHF. However, the TDLDA Rydberg excitation energies are much too low and pollute the valence spectrum. This is related to the fact that the LDA ionization energy, given by the opposite of the HOMO orbital (I = -ǫ HOMO ), is much too low. The somewhat opposite behaviors of TDHF and TDLDA suggest that more accurate valence and Rydberg excitation energies can be obtained with a proper combination of these two approaches. This is what is effectively done in TDDFT using (range-separated) hybrid functionals.

State

Transition TDHF TDHF-TDA TDLDA TDLDA-TDA EOM-CCSD Expt Valence excitation energies (eV)

κ

  ar κ * ar |Φ 0 + • • • ,

  ar,bs + B ar,bs ) ∂κ bs ∂E j = r|μ j |a .

n

  (t → -∞) = 0. Performing the integration over t ′ gives c(10)

  where A ar,bs = Φ r a | Ĥ0 -E HF |Φ s b = (ε r -ε a ) δ ab δ rs + rb||as , B ar,bs = Φ rs ab | Ĥ0 |Φ 0 = rs||ab and S ar,bs = Φ r a |Φ s b = δ ab δ rs .

  HF |Φ s b and B ar,bs = Φrs ab | Ĥ0 -E HF |Φ 0 = Φ rs ab | Ĥ0 |Φ 0 .

	2 E HF ∂κ ar ∂κ bs κ=0	= 2 (A ar,bs + B ar,bs ) ,	(36)
	with A ar,bs = Φ r a | Ĥ0 -E

Exercise: Using the Slater-Condon rules for calculating expectation values of operators over Slater determinants, show that

Table 1 :

 1

		UCHF CPHF Estimated exact
	He	1.00	1.32	1.383
	Ne	1.98	2.38	2.669
	Ar	10.1	10.8	11.08
	Kr	15.9	16.5	16.79
	Be	30.6	45.6	37.76
	Mg	55.2	81.6	71.26
	Ca	125	185	157.1

Spherically averaged static dipole polarizability ᾱ = (α x,x + α y,y + α z,z )/3 (in bohr 3 ) for rare-gas and alkaline-earth-metal atoms calculated at the UCHF and CPHF levels with uncontracted d-aug-cc-pCV5Z basis sets. Table extracted from Ref.

[START_REF] Toulouse | [END_REF]

.

Table 2 :

 2 Excitation energies and oscillator strengths of the CO molecule calculated by TDHF and TDLDA, without and with the Tamm-Dancoff approximation (TDA), and by EOM-CCSD, using the Sadlej+ basis set. Excitation energies above the ionization threshold are indicated in italics. Mean absolute deviations (MADs) are calculated with respect to the EOM-CCSD values. Experimental values are also given for comparison. Adapted from Ref.[6].

	3 Π	5σ → 2π *	5.28	5.85	5.95	6.04	6.45	6.32
	3 Σ +	1π → 2π *	6.33	7.79	8.38	8.54	8.42	8.51
	1 Π	5σ → 2π *	8.80	9.08	8.18	8.42	8.76	8.51
	3 ∆	1π → 2π *	7.87	8.74	9.16	9.20	9.39	9.36
	3 Σ -	1π → 2π *	9.37	9.73	9.84	9.84	9.97	9.88
	1 Σ -	1π → 2π *	9.37	9.73	9.84	9.84	10.19	9.88
	1 ∆	1π → 2π *	9.96	10.15	10.31	10.33	10.31	10.23
	3 Π	4σ → 2π *	13.05	13.31	11.40	11.43	12.49
				Rydberg excitation energies (eV)
	3 Σ +	5σ → 6σ	11.07	11.18	9.55	9.56	10.60	10.40
	1 Σ +	5σ → 6σ	12.23	12.27	9.93	9.95	11.15	10.78
	3 Σ +	5σ → 7σ	12.40	12.42	10.26	10.26	11.42	11.30
	1 Σ +	5σ → 7σ	12.78	12.79	10.47	10.50	11.64	11.40
	3 Π	5σ → 3π	12.52	12.60	10.39	10.39	11.66	11.55
	1 Π	5σ → 3π	12.87	12.88	10.48	10.50	11.84	11.53
				Ionization threshold: -ǫ HOMO (eV)
			15.11	15.11	9.12	9.12		15.58
		MAD of excitation energies with respect to EOM-CCSD (eV)
	Valence		0.89	0.49	0.37	0.33
	Rydberg		0.93	0.97	1.21	1.19
	Total		0.91	0.69	0.73	0.70
				Oscillator strengths (×10 -2 )
	1 Π	5σ → 2π *	8.55	11.48	8.69	11.41	8.66
	1 Σ +	5σ → 6σ	10.58	10.49	1.84	2.36	0.58
	1 Σ +	5σ → 7σ	9.39	10.22	12.53	12.25	20.71
	1 Π	5σ → 3π	5.13	4.94	2.71	2.16	4.94


Renaud Hedouin

Olivier Commowick

Julie Coloigner

Imaging biomarker discovery in major depressive disorder with diffusion MRI multi-compartment models

Introduction:

Major depressive disorder (MDD) is a disease widespread all over the world associated with a large and increasing economic, societal, and personal burden [Ferrari 2013]. Yet, the pathophysiology underlying MDD is not well understood, and rates of either non-response to treatment or MDD recurrence remain high. Many studies have already demonstrated impaired WM integrity, which corresponds to a FA decrease in the patient group compared to control subjects, with significant results found in several brain regions [Coloigner 2019]. Yet, the comparison between treatment-resistant and non-treatment-resistant populations has never been addressed. The FA at the voxel level provided by the classic diffusion tensor imaging (DTI) does not accurately describe the underlying microstructure. We propose a tract-based pipeline analysis using multi-compartment models (MCM) that allows to model water diffusion under the voxel resolution.

Methods:

We used the AMYNET database, constituted of 20 depressed unipolar resistant(TRD), 9 depressed unipolar in remission (nTRD) and 19 healthy subjects between 18 and 35 years. Half-noddi DWI sequences (similar The DTI model comports limitation, in particular, to disentangle complex sources of water diffusion, as crossing fiber. To tackle this issue, a MCM that decomposes the water diffusion as a weighted sum of 3 anisotropic tensors which represent fiber bundles in different directions along with 2 isotropic compartments for free water (FW), and iso-restricted water was estimated A deterministic tractography was performed on the entire brain atlas and filtered using TractSeg

[Wasserthal 2018] resulting in 72 identified fiber bundles. For all MCMs on the MNI atlas, each fiber bundle was used to extract microstructure properties for all patients using the closest anisotropic compartment (FA, MD, axial, radial) and the isotropic compartments (FW, iso-restricted) of the MCM. Finally, a two-sample t-test (Student's test) statistic analysis was performed along the tracts between the two groups.

Results:

We tested this entire pipeline on the AMYNET database. We compare in Figure 1, the FW weight difference between the TRD group (19 subjects) and the healthy group (19 subjects) along the left parieto-occipital pontine (POPT) tract. We can notice an increase of the FW within the TRD group compared with the healthy group with a significant p-value. This is consistent with previous studies linking degradation of WM, expected within the TRD population, with an increase of the FW weight [Jacobs 2015]. Let notice that those MCM parameters offer an interpretation of the brain structure at a microscopic level than cannot provide DTI parameters.

  [Stamm 2016]. To perform group comparison, all data were registered on the MNI152 atlas (an average of 152 T1-weighted MRI scans) [Grabner 2006] using dense transformation [Commowick 2012] estimated between the T1-weighted subject and T1-weighted atlas images. The transformations were directly applied to the MCMs using a specific registration framework [Hedouin 2020].

  Zhang 2012] while avoiding the redundancy of information) were acquired with 2 mm isotropic resolution, 97 gradients orientations distributed on 2 spheres (b = 700 and b = 3000 s.mm-1) and 7 b0. 3D MP2RAGE sequences were used to reconstruct a bias-free T1-weighted image for registration purposes [Marques 2010]. The data were processed with an open-source medical imaging toolbox Anima including eddy current correction, brain extraction, block-matching distortion correction [Hedouin 2017], NL-means denoising.

	1	1	1

1 OHBM https://ww3.aievolution.com/hbm2201/index.cfm?do=ab... 1 of 4 12/20/21, 15:12 to the noddi sequence [

Conclusions:

We proposed a method to treat high-quality data that can reveal subtle microstructure change using MCM. We showed an example of microstructure parameters tracts unavailable using the standard DTI. Further studies with larger populations including nTRD will be conducted. We will do so using MCM and quantitative MRI (relaxometry) in addition to clinical data using these multiple modalities as a whole.

Disorders of the


Table 1

Definitions and specific terms.

English terms Definition

Biomechanics Application of mechanical principles to living organisms Force Mechanical action resulting in the displacement, acceleration, or deformation of a solid Moment arm Distance from the axis of rotation to the applying point of the force Torque Ability of a force to rotate a solid Loads and constraints Forces and moments applied to the solid as external actions (Von Mises) stress Inner body pressure defined as the ratio between the applied force and the surface (MPa) Contact pressure Mechanical interaction due to contact forces between two solids (MPa) Boundary conditions Actions defined prior to the simulation, imposing a value on displacements or forces (ex: fixed tibial bone) Region of interest (ROI) Anatomical landmarks or surface on which the output data from the computational model are collected Among the mechanical conditions studied after TKA, a squatting movement is typically used in computer simulations [START_REF] Hoshi | Mobile-bearing insert used with total knee arthroplasty does not rotate on the tibial tray during a squatting activity: a cross-sectional study[END_REF][START_REF] Mizu-Uchi | Patient-specific computer model of dynamic squatting after total knee arthroplasty[END_REF][START_REF] Huang | Patellofemoral kinematics during deep knee flexion after total knee replacement: a computational simulation[END_REF]. The patient's ability to flex their knee as much as possible from an extended leg position is often considered representative of activities of daily living and correlated with clinical results [START_REF] Latypova | Patellar bone strain after total knee arthroplasty is correlated with bone mineral density and body mass index[END_REF]. An increase in patellar loads at about 60 • flexion has been observed in simulations during squatting after TKA [START_REF] Latypova | Patellar bone strain after total knee arthroplasty is correlated with bone mineral density and body mass index[END_REF][START_REF] Latypova | A patient-specific model of total knee arthroplasty to estimate patellar strain: a case study[END_REF]. However, the role of tension in the patellar retinaculum is unspecified in these simulations.

Various finite element models (FEM) have been described in the literature, each having their own complexity, precision and functionality [START_REF] Bori | Development and validation of an in-silico virtual testing rig for analyzing total knee arthroplasty performance during passive deep flexion: a feasibility study[END_REF][START_REF] Woiczinski | Development and validation of a weight-bearing finite element model for total knee replacement[END_REF][START_REF] Baldwin | Dynamic finite element knee simulation for evaluation of knee replacement mechanics[END_REF][START_REF] Innocenti | Contact forces in several TKA designs during squatting: a numerical sensitivity analysis[END_REF]. Some computer simulations have been proven effective at predicting the tibiofemoral loads during a squat [START_REF] Fitzpatrick | Influence of patellofemoral articular geometry and material on mechanics of the unresurfaced patella[END_REF], but few models have focused on evaluating the patellofemoral joint during flexion. The boundary conditions and loading conditions are very different in these various models, making them difficult to compare [START_REF] Viceconti | Extracting clinically relevant data from finite element simulations[END_REF][START_REF] Erdemir | Considerations for reporting finite element analysis studies in biomechanics[END_REF]. The aim of this study was to validate a new computational model driven by forces and moments and to analyze patellofemoral reaction forces and stress during squatting after TKA implantation (Table 1).

Materials and methods

Knee and implant geometry

Bone geometries were extracted from preoperative CT-scan slices (General Electric Healthcare, Milwaukee, WI, USA) from a single patient with end-stage knee osteoarthritis and varus deformity, operated on for a TKA (Fig. 1). Three-dimensional (3D) reconstructions were completed using open-source 3D image processing software (3D Slicer, version 4.11) based on manual segmentation of the femur, tibia, fibula and patella. The final parametric model was exported in STL format and included the volumes of each bone. The mechanical and anatomical axes and the lengths of the bone segments were collected from preoperative biplanar radiographs (EOS Imaging, Paris, France). Discretization of the model was done using fine mesh generation with tetrahedric elements that averaged 3 mm in size. This mesh was composed of 67,254 elements and 120,646 nodes (Fig. 2). The geometry of a posterior-stabilized TKA with fixed bearing was utilized using CAD software (SolidWorks, Dassault system, France).

Total knee arthroplasty implantation

The CAD software was then used to implant the TKA components (femoral component, tibial tray, tibial insert, polyethylene patellar dome component) using a mechanical alignment technique as follows: femoral and tibial bone cuts perpendicular to the femoral and tibial mechanical axes, respectively, 7 • of femoral valgus, 3 • of femoral external rotation, and 7 • of tibial slope. The implantation was performed by a surgeon specialized in total joint replacement to optimize the size and position of the components based on the TKA implanted in the same patient in vivo. A 10-mm fixed bearing polyethylene tibial insert was inserted between the femoral and tibial components. Implants were rigidly fixed and kinematically constrained to their bone, with no modelling of the cement fixation.

Material properties

The ligament structures were simulated as wire bodies aiming to reproduce physiological behavior (element working in tension but not in compression). The ligament attachments were defined using the initial imaging and validated by the surgeon (Fig. 2). The bones and implants were considered as rigid deformable bodies with fully linked interfaces (except for the patellofemoral and tibiofemoral joints). Friction coefficients of = 0.02 and = 0.05 were applied to the patellofemoral and tibiofemoral joints, respectively [START_REF] Huang | Stress distribution of the patellofemoral joint in the anatomic V-shape and curved dome-shape femoral component: a comparison of resurfaced and unresurfaced patellae[END_REF]. The mechanical properties attributed to materials (Table 2) and soft tissues (Table 3) were assumed to be homogeneous and isotropic, and derived from the literature.

Loading and boundary conditions

The Oxford Knee Rig setup [START_REF] Bourne | A functional analysis of various knee arthroplasties[END_REF], a well-controlled and reproducible experimental protocol, was reproduced in silico, while keeping the tibia fixed during flexion (Fig. 3). A constant 130 Newtons (N) load was applied to the center of the femoral head along an axis joining the centers of the femoral head and ankle, to simulate the loading conditions during a squatting motion while standing on both legs [START_REF] Huang | Stress distribution of the patellofemoral joint in the anatomic V-shape and curved dome-shape femoral component: a comparison of resurfaced and unresurfaced patellae[END_REF]. As such, the direction of force applied varied depending on the flexion angle, by inducing a flexion moment → C in the knee replacement (Fig. 3). The components Fx along

C = F × d d = L1 × sin ˇ1  = ˇ1 + ˇ2 d = L1 × sin  -ˇ2 Fx = F × cos ˇ2 and Fy = F × sin ˇ2
The parameters F, d, L1, ˇ1 and ˇ2 were defined as shown in Fig. 4. The quadriceps tendon was modeled using unidimensional viscoelastic elements (shocks and springs set in parallel). The damping coefficients were chosen to minimize dynamic effects. The action of the hamstring muscles was ignored during the flexion movement [START_REF] Bori | Development and validation of an in-silico virtual testing rig for analyzing total knee arthroplasty performance during passive deep flexion: a feasibility study[END_REF][START_REF] Victor | The influence of muscle load on tibiofemoral knee kinematics[END_REF]. The analysis was completed for each flexion position (20 • , 30 • , 40 • , 60 • , 80 • , 100 • ) after having first calculated the quadriceps stiffness, for which the knee was maintained in equilibrium due to external loads.

Validation study

The results were compared to experimental data (cadaver studies using the Oxford Knee Rig) or validated computer data with similar loading conditions. For a squatting motion, the comparability analysis used standardized forces of twice the load applied initially (noted BW=2*130N=260 N) as per Mason et al. [START_REF] Mason | Patellofemoral joint forces[END_REF]. The validation study included forces applied to the quadriceps tendon (F Q ) and patellar tendon (F PT ), contact forces (medially and laterally, and medial/lateral ratio) and contact pressure on the tibiofemoral joint, and the flexion angle at which post-cam contact occurred (defined by the flexion angle at which the contact force was no longer null during flexion).

Output data

Patellofemoral contact force (F PF , expressed in N) and Von Mises stress on the patellar implant (expressed in MPa) were evaluated for each position in flexion.

Results

Quadriceps force

Quadriceps force increased continuously with knee flexion angle, reaching a maximum of about 6 times the applied load at 100 • flexion (Fig. 5).

Patellar tendon force

The patellar tendon force F PT /quadriceps force F Q ratio was slightly reduced with increasing knee flexion angle (Fig. 6).

Tibiofemoral contact forces and pressure

Lateral and medial tibiofemoral contact forces increased continuously with knee flexion angle, up to 407 N medially and 257 N laterally at 100 • flexion. The medial-to-lateral distribution on the tibial insert was 64% and 36%, respectively. The contact pressure map on the tibial insert relative to knee flexion angle showed peak contact pressure more posteriorly with flexion, with a femoral rollback pattern (Fig. 7).

Post-cam contact forces

Post-cam contact forces started increasing slightly below 70 • of knee flexion.

Patellofemoral contact forces and Von Mises stress (output data)

Patellofemoral contact forces F PF increased from 20 • to 100 • , reaching 7 times the initial applied load (i.e. 3.5 times the BW value) (Fig. 8), and were correlated with increasing force applied to the quadriceps tendon (Fig. 9). The maximal F PF /F Q ratio was about 40% during our simulation. Von Mises stress on the patellar implant increased from 5.5 MPa (20 • flexion) to 16.4 MPa (100 • flexion) on average (Fig. 10). While they were initially centered, the stress shifted proximally and medially on the patellar implant, likely attributable to patellar contact with the intercondylar notch.

Discussion

The aim of this study was to validate a new finite element model of a TKA that provides reliable biomechanical data on the patellofemoral joint reaction forces during a simulated squatting motion. Finite element analysis is a method known for its precise calculation of joint loads, generating mechanical output data related to deformable elements [START_REF] Bori | Development and validation of an in-silico virtual testing rig for analyzing total knee arthroplasty performance during passive deep flexion: a feasibility study[END_REF]. Thus, this is a reference method for understanding the forces exerted on the extensor mechanism during flexion.

The initial step of the computer validation focuses on the tibiofemoral joint's loads and kinematics, similar to the in vivo or experimental conditions observed after implantation of posteriorstabilized TKA: presence of posterior femoral roll-back of several millimeters without paradoxical anterior translation, limited lateral lift-off of femoral condyles during flexion [START_REF] Argenson | In vivo kinematic evaluation and design considerations related to high flexion in total knee arthroplasty[END_REF][START_REF] Belvedere | In vivo kinematics of knee replacement during daily living activities: condylar and post-cam contact assessment by three-dimensional fluoroscopy and finite element analyses: contact assessment in total knee replacement[END_REF]. In our simulation, the tibiofemoral contact forces, their distribution between the medial and lateral compartments (about 60% vs. 40%), along with the instant of post-cam contact (typically reported to be between 60 • and 75 • flexion) were consistent with several published studies [START_REF] Argenson | In vivo kinematic evaluation and design considerations related to high flexion in total knee arthroplasty[END_REF][START_REF] Wang | The biomechanical effect of different posterior tibial slopes on the tibiofemoral joint after posterior-stabilized total knee arthroplasty[END_REF][START_REF] Innocenti | Contact forces in several TKA designs during squatting: a numerical sensitivity analysis[END_REF].

The innovative method of our approach included an a priori calculation of the quadriceps stiffness to equilibrate the model at the given flexion angle. The 130 N applied load was similar to the loading conditions used by Victor et al. [START_REF] Belvedere | In vivo kinematics of knee replacement during daily living activities: condylar and post-cam contact assessment by three-dimensional fluoroscopy and finite element analyses: contact assessment in total knee replacement[END_REF]. By standardizing external forces to twice the applied load, forces applied to the extensor mechanism were consistent with Victor et al. and other experimental studies in which external forces were determined by inverse kinematics [START_REF] Sharma | In vivo patellofemoral forces in high flexion total knee arthroplasty[END_REF][START_REF] Churchill | The influence of femoral rollback on patellofemoral contact loads in total knee arthroplasty[END_REF]. Quadriceps force was 6 times higher than the initially applied load at 100 • flexion, which is suitable for squatting motion [START_REF] Belvedere | In vivo kinematics of knee replacement during daily living activities: condylar and post-cam contact assessment by three-dimensional fluoroscopy and finite element analyses: contact assessment in total knee replacement[END_REF]. However, a large variability of quadriceps force values exists in the literature, ranging from one to two times our results [START_REF] Woiczinski | Development and validation of a weight-bearing finite element model for total knee replacement[END_REF]. Previously collected kinematic tibiofemoral data from experimental study can be incorporated and used to drive the finite element simulation, which can over-estimate quadriceps force [START_REF] Nakamura | Superiorinferior position of patellar component affects patellofemoral kinematics and contact forces in computer simulation[END_REF].

Patellofemoral contact force increased gradually during knee flexion and proportionally to the quadriceps force. Finally, the distribution pattern and magnitude of Von Mises stress on the patellar implant were comparable to the literature [START_REF] Sharma | In vivo patellofemoral forces in high flexion total knee arthroplasty[END_REF][START_REF] Churchill | The influence of femoral rollback on patellofemoral contact loads in total knee arthroplasty[END_REF][START_REF] Latypova | Patellar bone strain after total knee arthroplasty is correlated with bone mineral density and body mass index[END_REF]. Interestingly, the highest mean Von Mises stress found in our simulation ( This study has some limitations. First, a single implant design was used, thus these results cannot be generalized to other TKA implants. Second, the validation analysis was restricted to published studies with comparable loading and boundary conditions. While two types of mesh were tested with no major impact on the results, conducting a numerical sensitivity analysis by changing the mesh or material parameters could enhance the validity of our model. Additionally, our simulation was designed to ana- lyze extensor mechanism reaction forces and patellofemoral stress as output data, thus patellofemoral kinematics were not collected. Lastly, model inputs such as the patellar retinaculum [START_REF] Wang | Effect of patellar tendon release on squatting characteristics of tibiofemoral joint after total knee arthroplasty: a simulation analysis[END_REF][START_REF] Koh | Effect of post-cam design for normal knee joint kinematic, ligament, and quadriceps force in patient-specific posterior-stabilized total knee arthroplasty by using finite element analysis[END_REF], hamstring tendons [START_REF] Stevens-Lapsley | Quadriceps and hamstrings muscle dysfunction after total knee arthroplasty[END_REF] and more complex collateral ligaments models [START_REF] Beidokhti | The influence of ligament modelling strategies on the predictive capability of finite element models of the human knee joint[END_REF] could be incorporated to provide more realistic patellofemoral behavior during squatting motion. In addition to some intraoperative tools [START_REF] Ogawa | Extra-articular factors of the femur and tibia affecting knee balance in mechanically aligned total knee arthroplasty[END_REF], a preoperative computational simulation tool could predict individual patellofemoral forces and contact pressure with respect to implant position or patient-specific anatomy, providingdata that would help the surgeon in planning the TKA implantation.

Conclusion

This study showed agreement among litteratureand our computational model for patellofemoral forces and stress under squatting conditions after TKA implantation. This model furthers our understanding of the influence of implantation techniques on patellofemoral reaction forces during flexion.

Disclosure of interest

The authors declare that they have no competing interest.

Fig. 1 .

 1 Fig. 1. Methodology using CT-scan for bone reconstruction and implant modeling.

Fig. 2 .

 2 Fig. 2. Computational model of the total knee arthroplasty and the extensor mechanism using finite element analysis and mesh generation.

  → x0 and Fy along → y0 and the force → F and moment → C were determined using the following equations:

Fig. 3 .

 3 Fig. 3. Boundary conditions in the sagittal plane during a squat simulation. Initial load on the femoral head (A), expression of this load in the mobile femoral frame with flexion angle (B), definition of a new fixed knee frame (C), application of the initial load in the knee frame inducing a moment at the center of the knee replacement (D).

Fig. 4 .

 4 Fig. 4. Expression of the force

  16 MPa at 100 • flexion) was similar to experimental results from Glaser et al. (17 MPa at 90 • flexion), values surpassing the polyethylene yield point of the patellar implant (14.4 MPa) from 90 • of flexion [30].

Fig. 5 .

 5 Fig. 5. Changes in quadriceps tendon force relative to body weight (FQ/BW=) under squat loading conditions, compared to literature.

Fig. 6 .

 6 Fig. 6. Changes in patellar tendon force FPT / quadriceps force FQ ratio during flexion, compared to literature.

Fig. 7 .

 7 Fig. 7. Contact stress distribution on the tibial insert during flexion.

Fig. 8 .

 8 Fig. 8. Changes in patellofemoral contact forces relative to body weight (FPF/BW) under squat loading conditions, compared to literature. BW = twice the load applied initially (260 N)

Fig. 9 .

 9 Fig. 9. Changes in patellofemoral contact force FPF / quadriceps force FQ ratio , compared to literature.

Fig. 10 .

 10 Fig. 10. Von Mises stress distribution in the patellar implant during flexion.

Table 2

 2 TKA component properties.

	Structure	Material	Young's Modulus E (MPa)	Poisson coefficient
	Tibial insert, patellar implant	Polyethylene	685	0.4
	Femoral implant	CrCoMo	220,000	0.3
	Tibial implant	Titanium	110,000	0.3

Table 3

 3 Ligament properties.

	Structure	Young's modulus E (MPa)	Poisson coefficient	Surface (mm 2 )
	Lateral collateral ligament	180	0.4	9.42
	Medial collateral ligament	160	0.4	9.42
	Patellofemoral ligament	180	0.4	6.28
	Patellar tendon	1600	0.4	9.42

Funding

Agence Nationale de la Recherche (ANR): LabEx NuMEV Scholarship (SM).

Author contributions

L.D.: investigation, methodology, formal analysis, visualization, writing -original draft, writing -review and editing.

S.M.: investigation, methodology, computational analysis, visualization, writing-original draft, writing -review and editing. E.M.: writing -review and editing. F.C.: visualization, writing -review and editing. F.J.: conceptualization, computational analysis, project administration, writing -review and editing.