James Mure-Dubois

Heinz Hügli

Fusion of Time of Flight Camera Point Clouds

Recent time of flight cameras deliver range images (2.5D) in realtime, and can be considered as a significant improvement when compared to conventional (2D) cameras. However, the range map produced has only a limited extent, and suffers from occlusions. In this paper, we investigate fusion methods for partially overlapping range images, aiming to address the issues of lateral field of view extension (by combining depth images with parallel view axes) and occlusion removal (by imaging the same scene from different viewpoints).

Introduction

Recent time of flight (TOF) cameras enable new applications by producing range maps in real-time. Nevertheless, many situations arise where a single range map does not provide enough information. In this paper, we investigate real-time fusion of multiple range images acquired with TOF cameras. This fusion allows to obtain more information on the 3D scene observed. In particular, we will discuss the possibility to extend the field of view and to remove occlusions when using multi-camera networks.

In the following section, the operation principle of a TOF camera is described, and limitations of monocular view are presented. Then, available options to build TOF camera networks are discussed. Section 2 introduces two candidate calibration strategies for such networks : bundle adjustment and surface matching. Those methods are compared in section 3. Although the evaluation is only qualitative, the comparison clearly illustrates the advantages of the surface matching technique. Finally, section 4 presents example applications of multi-camera networks, compatible with real-time operation.

Operation of a Time of Flight Camera

Time of flight cameras involve active illumination, and deliver range (or depth) data by measuring the time needed for a light signal to travel from the camera light source to the scene and back to the camera sensor, as illustrated in figure 1. Present cameras ( [1], [2]) are based on the continuous emission of a periodic signal. The frequency of modulation f of this signal is typically 20MHz. The periodic signal S(i) received at each pixel (i) of the camera sensor is described by its amplitude A(i) and its phase ϕ(i). The range r is directly proportional to the phase. If we note c as the speed of light, we have :

S(i) = A(i) • e ϕ(i) . r(i) = c 4π f • ϕ(i) . (1)
One key advantage of such cameras is their real-time capability : range maps can be delivered at 20 frames per second (fps).

Single-Camera Operation, and its Limitations

Our test TOF device, the SwissRanger SR-3000 [1] has a 176 × 144 sensor. Figure 2 illustrates the range map r (i, j) produced by this device for a simple scene. By using a pinhole camera perspective projection model [START_REF] Tsai | A versatile Camera Calibration Technique for High-Accuracy 3D Machine Vision Metrology Using Off-the-Shelf TV Cameras and Lenses[END_REF], it is possible to transform this range map into a cloud of 3D points. We call f the camera focal length, d x and d y the pixel pitch in the x (resp. y) direction, and define normalized pixel coordinates (X c , Y c ), which are pixel coordinates relative to the position (c x , c y ) of the optical center on the sensor array. Neglecting lens distortion, the transformation between the range map r (i, j) and 3D coordinates (relative to the camera position) is given by :

z = r • f √ f 2 +(Xcdx) 2 +(Ycdy) 2 . x = z • Xcdx f . y = z • Ycdy f . (2)
The SR-3000 driver software includes functions performing this transformation from the range map r to the associated point cloud P, which can be rendered in a 3D visualization software, such as Paraview [START_REF]Kitware: Paraview -Parallel visualization application[END_REF] (fig. 3). This representation illustrates two limitations of range imaging with a TOF camera. First, the lateral field of view (FOV) is limited. In our example, only the top part of the person is included in the image. The second limitation comes from the perspective projection. Objects close to the camera occlude objects farther away. This effect is clearly illustrated by the shadow cast by the person over the wall in fig. 3. We will see in the following sections that multi-camera networks allow to overcome those limitations.

Multi-Camera Operation

Since TOF cameras are active devices with their own illumination source, special care must be taken when multiple cameras are operated simultaneously. In particular, it is desirable to avoid interference between different devices, which could lead to erroneous range measurements. In this paper, we will discuss only interference minimization for the SR-3000 camera [1], which provides range data based on the phase measurement of a periodic signal. Büttgen [START_REF] Büttgen | Extending time-of-flight optical 3D-imaging to extreme operating conditions[END_REF] gives a good overview of multi-camera operation options for this device. Options allowing to image the same scene with multiple cameras, without requiring an explicit synchronization of the different devices are Frequency Division Multiple Access (FDMA) and Code Division Multiple Access (CDMA). Although a method using CDMA has been proposed for SwissRanger cameras [START_REF] Büttgen | Extending time-of-flight optical 3D-imaging to extreme operating conditions[END_REF] , its implementation has not yet been released to camera users. Therefore, our experiments were carried out using the FDMA approach. Specifically, we operated a first camera at f 0 = 20 MHz and a second camera at f 1 = 21 MHz. Lange [START_REF] Lange | 3D Time-of-flight distance measurement with custom solid-state image sensors in CMOS/CCD-technology[END_REF] showed that under these conditions, the crosstalk between the two devices is less than 40 dB, as long as the integration time exceeds T int,min = 100 µs. Since SR-3000 cameras support 4 different operating frequencies, networks with 4 different SR devices can be used with the FDMA approach. For network with more cameras, the CDMA approach must be employed. Büttgen [START_REF] Büttgen | Extending time-of-flight optical 3D-imaging to extreme operating conditions[END_REF] shows that 15bit long pseudo-noise sequences can be used to define codes allowing to operate up to 1800 cameras without significant interference.

Calibration Strategy

In this section, we consider the problem of aligning two point clouds P 0 and P 1 acquired by TOF devices C 0 and C 1 , in the specific case where the two fields of view overlap partially. More specifically, the alignment procedure should allow to obtain the rigid body transformation T C0,C1 between the point of view of cameras C 0 and C 1 . This rigid body transformation can be decomposed in a rotation R and a translation T, and is therefore fully specified by 6 parameters : the three translation coefficients in T, and the rotation angles θ (yaw), φ (pitch) and ψ (tilt). If those parameters are known, it becomes possible, for each point x 1 in the cloud P 1 , to compute the corresponding position x 0 in the coordinates system of C 0

x 0 = T C0,C1 [x 1 ] = R • x 1 + T . (3)
We do not explicitly consider here the general case of global alignment for a network with more than 2 cameras. Nevertheless, the same approach can be used to align an arbitrary number of cameras, provided that a chain of pairwise overlapping FOV links all cameras in the network. Registration of multiple views of the same scene has been thoroughly investigated for standard (intensity, color) 2D cameras [START_REF] Tsai | A versatile Camera Calibration Technique for High-Accuracy 3D Machine Vision Metrology Using Off-the-Shelf TV Cameras and Lenses[END_REF][7] [START_REF] Bouguet | Camera calibration toolbox for matlab[END_REF], where very accurate results can be obtained by means of bundle adjustment. We discuss the applicability of this approach for point clouds produced with SR-3000 cameras. Then, we will present an ad-hoc registration procedure, based on the range images measured. This method is better suited to the SR-3000 specifications.

Bundle Adjustment

A classical approach based on a planar target was proposed by Tsai [START_REF] Tsai | A versatile Camera Calibration Technique for High-Accuracy 3D Machine Vision Metrology Using Off-the-Shelf TV Cameras and Lenses[END_REF] and later refined by Zhang [START_REF] Zhang | Flexible Camera Calibration By Viewing a Plane From Unknown Orientations[END_REF] and Bouguet [START_REF] Bouguet | Camera calibration toolbox for matlab[END_REF], at which point it was integrated into the OpenCV library [9]. In this approach, the planar target is imaged at different positions in the cameras' fields of view. Reference points (usually corners) are extracted in the images. Since the real-world geometry of the target object is known, the perspective projection model for the camera can be inferred from the acquired images by bundle adjustment for the feature points. This photogrammetric calibration procedure provides camera intrinsic parameters such as focal length, position of the optical center and distortion, but also extrinsic parameters, namely the position of the camera with respect to the target pattern origin. This approach can be used for the SR-3000 camera since an amplitude map A (i, j) is measured in addition to the range map r (i, j).

As noted above, this method allows to compute extrinsic parameters. If those parameters are matched for images of the same target pattern acquired by two devices C 0 , C 1 at different positions, the relative position of the two devices can be computed. It is important to note here that this alignment method is based only on intensity images produced by the TOF cameras : the TOF range map is not exploited in this alignment procedure. Kahlmann [START_REF] Kahlmann | Range imaging metrology : investigation, calibration and development[END_REF] notes that the sensor lateral resolution is too low to use standard calibration targets for bundle adjustment for the SR-3000 camera. Lindner et al. [START_REF] Lindner | Data Fusion and Edge-Enhanced Distance Refinement for 2D RGB and 3D Range Images[END_REF] successfully applied the OpenCV calibration procedure to register the range image of a PMD TOF camera with a color image produced by a standard camera (in this case, the lateral resolution of the color imager is higher than the TOF sensor resolution). We carried out experiments based on the calibration toolbox provided by Bouguet [START_REF] Bouguet | Camera calibration toolbox for matlab[END_REF], using a large checkerboard pattern (each square was 100 mm). Figure 4 shows an example of calibration images, along with a rendering of the computed relative camera positions. Section 3.1 illustrates that the extrinsic parameters obtained do not allow to successfully align point clouds obtained with two SR-3000 cameras. One of the key issue is that the camera model parameters computed from the amplitude images may be different from the parameters used in the SR-3000 driver software for the transformation of the range map into a 3D cloud of points. Therefore, we investigated a more robust method to align range images, described in the following section. This method is based on a simple calibration scene, and relies only on the camera model used for point cloud generation in the SR-3000 driver software.

Surface Matching Techniques

The method presented above is unsatisfactory since it requires the determination of a rigid body transformation between two camera viewpoints. Those viewpoints are usually far away from the measured point clouds, so that any error in the transformation parameters is greatly amplified. Moreover, only a small set of matched points are used to compute the transformation. A significant improvement can be achieved if most of the measured points are considered in the matching procedure. This is the case in automatic surface registration methods such as iterative closest points (ICP) [START_REF] Jost | A Multi-Resolution ICP with Heuristic Closest Point Search for Fast and Robust 3D Registration of Range Images[END_REF]. Unfortunately, the TOF data obtained from real-time cameras is very noisy. Construction of valid surface sets from the range map is complicated, since pixel noise generates large surface patches, which would prevent an automatic surface registration procedure from converging. Therefore, we did not further investigate those methods in the present work. Nevertheless, as noise levels and denoising techniques for TOF device will probably improve in the next few years, this solution should be carefully investigated in future works.

We propose here an hybrid method, more suited to the SR-3000 data. We explicitly separate the surface matching in two steps : matching of rotation, and matching of translation between the two point clouds. In order to match the rotation, a planar target is imaged simultaneously with both cameras. Then, a random sample consensus (RANSAC) [START_REF] Torr | The development and comparison of robust methods for estimating the fundamental matrix[END_REF] method is used to recognize this planar regions in both point clouds and to compute the normal n to this object in the coordinates system of each camera. This feature is relatively robust to measurement noise, since a large number of points are taken into account for its determination. Point clouds for each camera are subsequently rotated so that this normal becomes parallel to the z axis. Performing the rotation for both cameras ensures that calibration errors are evenly distributed between both clouds. When this step is completed, the remaining translation parameters are computed from a set of 6 to 10 point pairs defined by a human operator. Again, using more points helps in lowering the impact of measurement noise.

Calibration Experiments

Calibration experiments were carried out using an ad-hoc TOF image acquisition application, allowing to record range map simultaneously with two SR-3000 devices, operated at frequencies f 0 = 20 MHz and f 1 = 21 MHz respectively. In the experiments presented here, the aim was set to determining the coordinates transformation T C0,C1 allowing the computation of the position of points P 1 measured by the camera C 1 into the coordinates system of the reference camera C 0 . A qualitative evaluation of the calibration performance is therefore possible : point clouds P 0 and P 1 are rendered in the same 3D scene, allowing an human observer to evaluate the fitting of both datasets.

Camera Calibration Toolbox

This toolbox uses a planar checkerboard pattern as calibration target (fig. 4a,b), and is based solely on amplitude data. Prior to extrinsic parameters determination, this toolbox allows to determine intrinsic parameters of the camera, such as the focal length f or the coordinates (c x , c y ) of the principal point. In tables 1 and 2 we compare the values obtained with this toolbox to the manufacturer's data for two SR-3X00 devices. In both cases, the agreement is good : the manufacturer's data lies within the uncertainty domains of calibration values. Unfortunately, the extrinsic parameters obtained with this calibration method proved useless for point cloud registration (fig. 5). Although disappointing, those bad registration results were expected since the camera model for image registration is (slightly) different from the camera used in point cloud computation. Moreover, variations related to fixed pattern noise at different frequencies further degrade the fitting between the two point clouds.

RANSAC based Method

To overcome errors due to inaccurate rotation parameters, we used the alignment procedure proposed in section 2.2. For the RANSAC procedure, the number of iterations N was set to 1600 and the threshold for inliers was d = 650 mm. Although this threshold is quite large, the obtained rotation parameters do not suffer greatly from noise, since the plane orientation is determined by averaging on a sufficiently large number of inliers. Typically, good results are obtained when the number of inliers exceeds 5000.

The remaining translation degrees of freedom were handled by manually defining 10 point pairs in the amplitude images provided by the SwissRanger cameras. Successful point cloud fusion using this calibration method is illustrated in fig. 8. This example is discussed in the next section.

Comparison of Calibration Results

The efficiency of the two alignment method discussed above can be compared by employing both methods for the same scene. Table 3 provides a comparison of the rotation matrix R and translation vector T for the bundle adjustment method and the RANSAC based method. Large discrepancies can be observed between the two results. Comparison of fused point clouds allows to determine which method performs best. Figure 6 shows an example scene acquired simultaneously with two TOF cameras. As stated above, results obtained with camera calibration toolbox are clearly too inaccurate for successful point cloud merging (fig. 7). On the contrary, the hybrid method where rotation is determined by matching plane primitives obtained through RANSAC allows to successfully fuse the point clouds (fig. 8) : all objects in the scene are uniquely represented in the merged point cloud. In order to go beyond the qualitative discussion presented here, an appropriate distance function between two TOF point clouds must be defined. Work on this topic is currently in progress, and should ultimately allow to determine a least error transformation (R, T). We provide below preliminary quantitative results for the estimation of the accuracy of the calibration based on RANSAC planes. Accuracy of RANSAC based Calibration Since no ground truth data is available, the accuracy estimation must be carried out by quantifying the correspondence between the 3D point clouds merged. Laveoué et al. [START_REF] Lavoué | Perceptually driven 3D distance metrics with application to watermarking[END_REF] proposed a 3D quality metric related to human perception, but this metric is defined only for meshes. Work is currently in progress to define a quality metric adapted to data delivered by TOF camera networks.

Table 3. Comparison of coordinates transformation determined by different alignment methods

For preliminary results, we use position measurements of a simple target identified by its amplitude pattern. The target position was measured at 6 different stations. An example set of images for this procedure is given in figure 9. For each station, the distance from the target to the closest camera is larger than 2m. Table 4 summarizes the measured discrepancies for the target position after calibration. The average error is less than 60 mm, for all directions, but the variation is high, indicating that more tests will be required to confirm this result. Moreover, the observations are not centered around zero, indicating that a systematic error is still present. Nevertheless, the performance achieved is good, especially when compared to the typical noise level for a TOF camera: 150 mm at this range. Those preliminary results confirm that the RANSAC based calibration method performs well in a simple network. Future experiments will help to determine if systematic errors can be further reduced, and if the global accuracy is reduced when more cameras are present in the network.

Synchronous Imaging with Calibrated Cameras

The two problems we are aiming to address with synchronous multi-camera acquisition are extension of the field of view, and removal of occlusion. A TOF camera network allows to perform synchronous acquisition of sequences of range images for scenes containing motion. Parameters for point cloud registration are computed offline, through an interactive procedure using an initial transform provided by the RANSAC registration procedure described in sec. 3.2. The computed point cloud are then visualized in real-time in the same 3D renderer, allowing to produce movies of 3D scenes.

FOV Extension

Figure 10 illustrates field of view extension. In this experiment, a person walks in front of two TOF cameras. The reference camera C 0 is aimed at the upper body, while the second camera image toward the legs. The resulting point clouds clearly shows the reconstruction of the full person walking. We emphasize here that the fusion is a realtime operation. The example of figure 10 is just a frozen frame of a sequence recorded at 20 fps.

Occlusion Removal

Removal of occlusions is illustrated in fig. 11. Using a single TOF camera results in large unmonitored areas on the wall behind a person close to the camera (fig. 11b). The occlusion is removed when a second TOF camera is added (fig. 11d), allowing to confirm that a single human is present in the cameras field of view. A typical application of occlusion resolution is human monitoring applications, where the number of humans present in a defined area must be reliably determined.

Conclusions

Applications involving synchronous operation TOF cameras in a multi-camera network were investigated. Experiments carried out with SR-3000 cameras show that registration methods used for conventional cameras can not be applied directly to TOF imagers with coarse sampling. An hybrid method using RANSAC plane primitives for alignment was introduced. This method provides better results, and allows to fuse points clouds acquired by the TOF network. While calibration is performed offline, fusion is carried out in real-time, allowing the network to operate at the TOF camera native speed (20 fps). Field of view extension and occlusion removal were realized with multi-camera networks, for scenes with motion. Possible improvements to the present work include definition of an error metric for point cloud fusion allowing to determine a least error solution, and extensions to support this error minimization for an arbitrary number of TOF sensors in the network.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Time of flight camera -Principle of operation

Fig. 3 .

 3 Fig. 3. Cloud of points P obtained from the TOF range map in figure 2b.

Fig. 4 .

 4 Fig. 4. Bundle adjustment calibration toolbox results

( a )Fig. 5 .

 a5 Fig. 5. Point cloud combination with parameters computed by camera calibration toolbox. Large discrepancies between the point clouds can be observed.

Fig. 6 .

 6 Fig. 6. Amplitude and range images acquired simultaneously with two TOF devices.

  Left point cloud -Perspective view (b) Combined point cloud -Perspective view (c) Right point cloud -Perspective view (d) Left point cloud -Top view (e) Combined point cloud -Top view (f ) Right point cloud -Top view

Fig. 7 .

 7 Fig. 7. Point cloud fusion with bundle adjustment calibration. Large discrepancies between the two point clouds can still be observed.

Fig. 8 .

 8 Fig. 8. Point cloud fusion with RANSAC calibration. The fused point cloud provides more information on the scene than individual clouds.

( a )

 a Left camera -Amplitude (b) Right camera -Amplitude (c) Combined point clouds

Fig. 9 .

 9 Fig. 9. Error estimation by measuring the position of a simple target.

Fig. 10 .

 10 Fig. 10. Point clouds combination for field of view extension.

Fig. 11 .

 11 Fig. 11. Point clouds combination for occlusions removal.

Table 1 .

 1 Intrinsic camera parameters for SR-3100 (sn097027)

	Parameter MESA	calib.
	f [mm]	8.0 ±n.a. 8.04 ±0.23
	cx	85.0 ±n.a. 83.5 ± 4.7
	cy	76.7 ±n.a. 80.3 ± 5.4

Table 2 .

 2 Intrinsic camera parameters for SR-3000 (sn296012)

	Parameter MESA	calib.
	f [mm]	8.0 ±n.a. 7.98 ±0.26
	cx	95.1 ±n.a. 93.8 ± 5.2
	cy	56.3 ±n.a. 51.6 ± 5.9

Table 4 .

 4 Discrepancies in position of simple target (6 stations).

	Station	1	2	3	4	5	6	Avg std
	∆x [mm] +46.3 +70.2 +32.8 +56.7 +49.0 -09.6 +40.9 27.6
	∆y [mm] +65.2 +60.3 +61.9 +76.4 +39.0 +32.1 +55.8 16.8
	∆z [mm] +13.9 +02.1 +18.1 +20.4 +40.8 +64.4 +26.6 22.4


H Le Quang

Q.-C He

Compact explicit matrix representations of the flexoelectric tensor and a graphic method for identifying all of its rotation and reflection symmetries

Keywords:

Flexoelectricity is an electromechanical phenomenon produced in a dielectric material, with or without centrosymmetric microstructure, undergoing a non-uniform strain. It is characterized by the fourth-order flexoelectric tensor which links the electric polarization vector with the gradient of the second-order strain tensor. Our previous work [Le Quang and He, Roy. Soc. London, Ser. A 467, 2369-2386, 2011] solved the fundamental theoretical problem of determining the number and types of all rotational symmetries that the flexoelectric tensor can exhibit. In the present one, compact explicit matrix representations of the flexoelectric tensor are provided so as to facilitate the use of it with any possible rotational symmetry. The number and types of all reflection symmetries that the flexoelectric tensor can have are also determined. To identify the rotational symmetry and reflection symmetry of a given flexoelectric tensor, a simple and efficient graphic method based on the concept of pole figures is presented and illustrated.

I. INTRODUCTION

Flexoelectricity is a coupled electromechanical phenomenon appearing a dielectric material subjected to a non-uniform strain. In contrast to piezoelectricity, it can be generated even in a dielectric material whose microstructure is centrosymmetric. Indeed, in the case of a dielectric material undergoing a uniform strain, the electric polarization is produced if and only if the microstructure of this material is non-centrosymmetric. However, when a dielectric material with (or without) a centrosymmetric microstructure is subjected to a nonuniform strain whose gradient is non null, a relative displacement of the centers of the positive and negative charges is resulted in and gives rise to an electric polarization.

The flexoelectric effects can be produced in a multitude of situations, for example, in bending crystal plates [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] , nanobeams and nanowires [START_REF] Liu | Effect of flexoelectricity on electrostatic potential in a bent piezoelectric nanowire[END_REF][START_REF] Nguyen | Nanoscale flexoelectricity[END_REF][START_REF] Yan | Flexoelectric effect on the electroelastic responses of bending piezoelectric nanobeams[END_REF] or when stretching thin films [START_REF] Catalan | The effect of flexoelectricity on the dielectric properties of inhomogeneously strained ferroelectric thin films[END_REF] on liquid crystals [START_REF] Meyer | Piezoelectric effects in liquid crystals[END_REF] and on elastomers [START_REF] Marvan | Flexoelectric effect in elastomers[END_REF] . The flexoelectric constants describing the flexoelectric effects of some dielectric materials were observed and measured in a direct or indirect way in a few experimental works such as those made by Ma and Cross 8-10 and Zubko et al. [START_REF] Zubko | Strain-gradient-induced polarization in SrTiO 3 single crystals[END_REF] for various perovskites which exhibit unusually high flexoelectricity, the ones of Kalinin and Meunier [START_REF] Kalinin | Electronic flexoelectricity in low-dimensional systems[END_REF] and Naumov et al. [START_REF] Naumov | Unusual flexoelectric effect in two-dimensional noncentrosymmetric sp 2bonded crystals[END_REF] for low-dimensional structures like nanographitic systems and two-dimensional boronnitride sheets or by Zhang et al. [START_REF] Zhang | Shear flexoelectric response along 3121 direction in polyvinylidene fluoride[END_REF][START_REF] Zhang | Investigation of the 2312 flexoelectric coefficient component of polyvinylidene fluoride: Deduction, simulation, and mensuration[END_REF] , Chu and Salem [START_REF] Chu | Flexoelectricity in several thermoplastic and thermosetting polymers[END_REF] and Zhou et al. [START_REF] Zhou | Flexoelectric effect in PVDF-based polymers[END_REF] for dielectric materials and polymers a) Corresponding author Tel: 33 (0) 160 957 797; Fax: 33 (0) 160 957 799; Email: hung.le-quang@univ-eiffel.fr such as TiO 2 ceramics and the polyvinylidene fluoride (PVDF). In parallel with these experimental works, theoretical studies were also conducted to demonstrate the size-dependent flexoelectric properties and surface effect of dielectric materials/structures in nanoscale, for example, by Sahin and Dost 18 , Tagantsev [START_REF] Tagantsev | Piezoelectricity and flexoelectricity in crystalline dielectrics[END_REF][START_REF] Tagantsev | Electric polarization in crystals and its response to thermal and elastic perturbation[END_REF] , Yurkov and Tagantsev [START_REF] Yurkov | Strong surface effect on direct bulk flexoelectric response in solids[END_REF] , He et al. [START_REF] He | Size dependent flexoelectric and mechanical properties of barium titanate nanobelt: A molecular dynamics study[END_REF] , Qi et al. [START_REF] Qi | Size-dependent bending of an electro-elastic bilayer nanobeam due to flexoelectricity and strain gradient elastic effect[END_REF] and Bai et al. [START_REF] Bai | Size dependent flexocaloric effect of paraelectric Ba 0.67 Sr 0.33 TiO 3 nanostructures[END_REF] . Numerical approaches were elaborated according either to the first-principles method, for example, by Maranganti and Sharma [START_REF] Maranganti | Atomistic determination of flexoelectric properties of crystalline dielectrics[END_REF] , Hong et al. [START_REF] Hong | First-principles theory of frozen-ion flexoelectricity[END_REF][START_REF] Hong | First-principles theory and calculation of flexoelectricity[END_REF] or to the other theoretical calculation methods such as finite element method by Deng et al. [START_REF] Deng | A three-dimensional mixed finite element for flexoelectricity[END_REF] and Yvonnet et al. [START_REF] Yvonnet | Apparent flexoelectricity due to heterogeneous piezoelectricity[END_REF] , phase-field method by Li et al. [START_REF] Li | Quantification of flexoelectricity in PbTiO 3 /SrTiO 3 superlattice polar vortices using machine learning and phase-field modeling[END_REF] and Wang et al. [START_REF] Wang | Phase-field modeling and electronic structural analysis of flexoelectric effect at 180 o domain walls in ferroelectric PbTiO 3[END_REF] to estimate the flexoelectric properties of some dielectric materials/structures. For more references about flexoelectricity and for discussions on potential important applications of flexoelectricity, the reader is referred to Tagantsev et al. [START_REF] Tagantsev | Electric polarization in crystals and its response to thermal and elastic perturbation[END_REF] , Sharma et al. [START_REF] Sharma | On the possibility of piezoelectric nanocomposites without using piezoelectric materials[END_REF] , Zubko et al. [START_REF] Zubko | Flexoelectric effect in solids[END_REF] , Wang et al. [START_REF] Wang | Flexoelectricity in solids: Progress, challenges, and perspectives[END_REF] , Narvaez et al. [START_REF] Narvaez | Enhanced flexoelectric-like response in oxide semiconductors[END_REF] , Abdollahi et al. [START_REF] Abdollahi | Piezoelectric Mimicry of Flexoelectricity[END_REF] and Shu et al. [START_REF] Shu | Flexoelectric materials and their related applications: A focused review[END_REF][START_REF] Shu | Photoflexoelectric effect in halide perovskites[END_REF] .

When a dielectric material is subjected to small deformations and when the piezoelectric and flexoelectric phenomena produced in it are linear, the electric polarization vector p is related to the infinitesimal strain tensor ε and the gradient of the latter, namely E = ∇ε, by a linear relation:

p i = D ijk ε jk + F ijkl E jkl . (1)
Above, D ijk are the matrix components of the thirdorder piezoelectric tensor D and F ijkl stand for the matrix components of the fourth-order flexoelectric tensor F. Due to the symmetry ε ij = ε ji of ε, the strain gradient E possesses the property E ijk = E jik and the matrix components of F have the following index permutation symmetry:

F ijkl = F ikjl . (2)
Note that, if the microstructure of the dielectric material in question exhibits centrosymmetry, the requirement that the third-order tensor D be invariant under the central inversion transformation implies that D is null, so that the constitutive law (1) reduces to

p i = F ijkl E jkl . (3)
In other words, when the dielectric material has a centrosymmetric microstructure, the piezoelectric effect disappears and the electric polarization vector p depends only on the strain gradient E.

In the linear constitutive law (1), the classical thirdorder piezoelectric tensor D has been completely investigated and understood; however, the fourth-order flexoelectric tensor F, which is much more complicated than the usual fourth-order elastic tensor, is far from being thoroughly studied and understood. In our previous one [START_REF] Quang | The number and types of all possible rotational symmetries for flexoelectric tensors[END_REF] , the number and types of all possible rotational symmetries for the flexoelectric tensor F were specified and the number of independent material parameters of F belonging to each possible symmetry class was determined. Later, Shu et al. [START_REF] Shu | Symmetry of flexoelectric coefficients in crystalline medium[END_REF] gave the matrix representations of F for various symmetries.

The present work can be regarded as a continuation of our previous one [START_REF] Quang | The number and types of all possible rotational symmetries for flexoelectric tensors[END_REF] . Precisely, novel 3 × 18 matrix representations of the flexoelectric tensor F are provided for all 12 rotational symmetries determined in our previous study [START_REF] Quang | The number and types of all possible rotational symmetries for flexoelectric tensors[END_REF] . The matrix representations of F given in the present work are well-structured and much more compact than those of Shu et al. [START_REF] Shu | Symmetry of flexoelectric coefficients in crystalline medium[END_REF] . This should facilitate the practical use of F in various anisotropic cases. In addition, the flexoelectric tensor F is further investigated in the present work by determining all of its reflection symmetries. It is proved that the 12 rotational symmetry classes of the flexoelectric tensor F are reduced to 8 reflection symmetries and these 8 reflection symmetries are identical to those of the fourth-order elastic tensor. Finally, a simple but efficient graphic method based on the notion of pole figures is suggested and illustrated for identifying the reflection symmetry and rotational symmetry that a given flexoelectric tensor may have.

The paper is structured as follows. In Section II, some notations and definitions used throughout the paper are presented. In Section III, the main results obtained in our previous work [START_REF] Quang | The number and types of all possible rotational symmetries for flexoelectric tensors[END_REF] on the symmetry groups and symmetry classes of the flexoelectric tensor are recalled for the paper to be self-contained. In Section IV, compact explicit matrix representations of the flexoelectric tensor for all possible rotational symmetries are provided. Section V is dedicated to finding out the reflection symmetry classes of the flexoelectric tensor. In Section VI, a simple and efficient graphic method is elaborated to identify not only the reflection symmetry but also the rotational symmetry of a given flexoelectric tensor. In Section VII, some concluding remarks are given.

II. NOTATIONS AND DEFINITIONS

Let V be a three-dimensional (3D) inner-product space over the reals R and Lin be the space of all linear transformations (second-order tensors) on V. The inner product of two vectors a and b of V is symbolized by a • b. The 3D orthogonal group O(3) is defined as

O(3) = {Q ∈ Lin | Qa • Qb = a • b, ∀a, b ∈ V}. The 3D rotation group SO(3) is given by SO(3) = {Q ∈ O(3) | det Q = 1}.
In what follows, Q(a, θ) stands for the rotation about a ∈ V through an angle θ ∈ [0, 2π). In particular, Q, Q and Q denote, respectively, the rotations Q (e 1 + e 2 + e 3 , 2π/3), Q 2e 2 + ϕ 2 e 3 , 2π/3 and Q (ϕe 2 + e 3 , π) with the golden ratio ϕ = (1 + √ 5)/2. For later use, it is convenient to introduce the following standard group notations:

(i) the identity group, denoted by I, is formed by the second-order identity tensor I;

(ii) the cyclic group Z r (r ≥ 2) contains r elements generated by Q(e 3 , 2π/r). In particular, when r → ∞, Z r becomes the group SO(2) consisting of all rotations Q about e 3 such that Qe 3 = e 3 ;

(iii) the dihedral group D r (r ≥ 2) comprises 2r elements generated by Q(e 3 , 2π/r) and Q(e 1 , π). The corresponding form of D r when r → ∞ is O(2) consisting of all orthogonal tensors Q such that Qe 3 = ±e 3 ;

(iv) the spatial groups T , O and I with T representing the tetrahedral group of 12 elements generated by D 2 and Q, O being the octahedral group of 24 elements generated by D 4 and Q, and I symbolizing the dodecahedral group of 60 elements generated by D 5 , Q and Q. Recall that the spatial groups T , O and I map a tetrahedron, a cube and a dodecahedron onto themselves, respectively.

Next, we define the space of flexoelectricity tensors as follows:

F = {F = F ijkl e i ⊗ e j ⊗ e k ⊗ e l | F ijkl = F ikjl }.
The symmetry group of a flexoelectricity tensor F ∈ F, is denoted by G(F) and characterized as

G(F) = {Q ∈ SO(3) | Q * F = F} (4)
where

Q * F = Q ir Q jm Q kn Q ls F rmns e i ⊗e j ⊗e k ⊗e l . The definition of G(F) implies that G(F) is a closed subgroup of SO(3) and G(Q * F) = QG(F)Q T for any orthogonal tensor Q ∈O(3). On the other hand, F is said to exhibit G-symmetry when G ⊆ G(F) and Q * F = F for all Q ∈ G.
Consequently, two materials characterized by their respective flexoelectric tensors F 1 ∈ F and F 2 ∈ F are said to have the same type of symmetry if and only if the symmetry groups of F 1 and F 2 are conjugate to each other, i.e.

F 1 ∼ F 2 ⇔ G(F 1 ) ∼ G(F 2 ) ⇔ ∃Q ∈ SO(3) such that G(F 1 ) = QG(F 2 )Q T . (5)
(e3, π), Q(e1, π), Q(v, 2π/3) with v = √ 3 3 (e1 + e2 + e3) O
Octahedral group of 24 elements generated by

Q(e3, π/2), Q(e1, π), Q(v, 2π/3) with v = √ 3 3 (e1 + e2 + e3) I Dodecahedron group of 60 elements generated by Q(e3, 2π/5), Q(e1, π), Q(w, 2π/3) with w = 1 √ 4+φ 2 (2e2 + φe3) and φ = √ 5+1 2 SO(2) Subgroup of rotations Q(e3, θ) with θ ∈ [0; 2π) O(2)
Subgroup generated by SO(2) and Q(e1, π) P(n)

Reflection transformation through the plane of normal n P F Set of reflection symmetry elements of F {P F } Reflection symmetry class of F P h

Set containing one reflection P(e3)

Pv k (k ≥ 1) Set containing k elements Pv k = {P(r3( 2pπ k )} 1≤p≤k with r3(θ) = sin θe1 + cos θe2 P hv k (k ≥ 1)
Set containing k elements of Pv k completed by P(e3) PO Cubic set consisting of 9 reflections with respect to the nine planes of which the normals of 6 pass through the center of each edge of a regular cube and the normals of 3 through the center of each face of the latter PI Icosahedral set of 15 reflections with respect to the fifteen planes whose normals pass through the center of each edge of a regular icosahedron

P O(3)
Set composed of all reflections P(n)

With the above notion of conjugacy, a family of nonempty subsets, (F i ) 1≤i≤N , of F acts as a partition of the flexoelectric tensor space F in the sense that no two elements of (F i ) 1≤i≤N overlap and the union of (F i ) 1≤i≤N is equal to F. Thus, each element F i of this partition characterizes a symmetry class for flexoelectric tensors.

On the other hand, for a given flexoelectric tensor F i ∈ F i with the symmetry group G(F i ), the collection of all the conjugates of G(F i ) in the set of subgroups of SO(3), i.e.

{G i } = {G(F i )} = { G ⊆ SO(3) | G = QG(F i )Q T , Q ∈SO(3)}, (6)
constitutes an intrinsic characterization of the type of rotational symmetries exhibited by the elements of F i .

Clearly, the definition of F i through {G i } is more convenient. Finally, we denote by {G} the collection of all the conjugates of G ∈ SO(3) in the set of subgroups of SO(3) and define F(G) as the set

F(G) = {F ∈ F | G(F) ∈ {G}}. (7)
For the convenience of the reader, the notations used in this paper is summarized in Table I.

III. ROTATIONAL SYMMETRY CLASSES OF THE FLEXOELECTRIC TENSOR

For the paper to be self-contained and for later use, the present section recalls the main results of our previous study [START_REF] Quang | The number and types of all possible rotational symmetries for flexoelectric tensors[END_REF] concerning the determination of the number and types of all the rotational symmetries for the flexoelectric tensor. First, using a general method due to Spencer [START_REF] Spencer | A note on the decomposition of tensors into traceless symmetric tensors[END_REF] , any fourth-order flexoelectric tensor F ∈ F can be first decomposed into totally symmetric tensors and then split into harmonic tensors. The following explicit harmonic decomposition is established for the flexoelectric tensor F: [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] klm + jkm [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] ilm + ikm [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] jlm ) + 1 4

F ijkl = [H] ijkl + 1 3 (3 ijm [H (
( ilm [H [START_REF] Liu | Effect of flexoelectricity on electrostatic potential in a bent piezoelectric nanowire[END_REF] ] jkm + jlm [H [START_REF] Liu | Effect of flexoelectricity on electrostatic potential in a bent piezoelectric nanowire[END_REF] ] ikm + klm [H [START_REF] Liu | Effect of flexoelectricity on electrostatic potential in a bent piezoelectric nanowire[END_REF] ] ijm ) [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] kl + δ ik [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] jl + δ il [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] jk + δ jk [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] il + δ jl [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] ik + δ kl [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] ij )

+ 1 7 (δ ij [H (
+ 1 6 (3 ijm kln + ikm jln + jkm iln )[H (2) ] mn + 2 9 (2δ jk [H (3) ] il -δ ik [H (3) ] jl + δ jl [H (3) ] ik -2δ il [H (3) ] jk + δ lk [H (3) ] ji -δ ij [H (3) ] kl ) + 1 6 (δ jl [H (4) ] ik + δ kl [H (4) ] ij + δ il [H (4) ] jk -δ ij [H (4) ] kl -δ ik [H (4) ] jl -δ jk [H (4) ] il ) + 3 10 (δ ij klm + δ jk ilm + δ ik jlm )[a (1) ] m + 1 12 (2 kli [a (2) ] j -klj [a (2) ] i ) + 1 12 (2 kim δ jl + 2 lim δ jk -ljm δ ik -kjm δ il -lkm δ ij )[a (2) ] m + 1 15 (11 ikm δ jl -4 jkm δ il -5 jlm δ ik + 10 ilm δ jk -5 klm δ ij + 3 ijm δ kl )[a (3) ] m + 1 15 (3 ijk [a (3) ] l + 3 ijl [a (3) ] k + ikl [a (3) ] j + jkl [a (3) ] i ) + α 1 15 (δ ij δ kl + δ ik δ jl + δ il δ jk ) + α 2 3 (δ jl δ ik + δ ij δ kl -2δ il δ jk ). ( 8
)
It can be seen from ( 8) that the harmonic decomposition of F contains : two scalars α 1 and α 2 ; three vectors a [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] , a [START_REF] Liu | Effect of flexoelectricity on electrostatic potential in a bent piezoelectric nanowire[END_REF] and a [START_REF] Nguyen | Nanoscale flexoelectricity[END_REF] ; four second-order harmonic tensors H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] , H [START_REF] Liu | Effect of flexoelectricity on electrostatic potential in a bent piezoelectric nanowire[END_REF] , H [START_REF] Nguyen | Nanoscale flexoelectricity[END_REF] and H [START_REF] Yan | Flexoelectric effect on the electroelastic responses of bending piezoelectric nanobeams[END_REF] ; two third-order harmonic tensors H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] and H [START_REF] Liu | Effect of flexoelectricity on electrostatic potential in a bent piezoelectric nanowire[END_REF] and a fourth-order harmonic tensor H. The components of these harmonic tensors are explicitly expressed in terms of F ijkl as

α 1 = 1 3 (F pqqp + 2F ppqq ), α 2 = 1 3 pqk mnk F (mp)nq , (9)
[a [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF]

] k = 1 9 pqk (F pmmq + 2F mmpq ), [a (2) ] k = mnp F mnkp , [a (3) ] k = 1 6 ( pqk F pqmm + pqn F pqkn ), (10)
[H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] km = F (ppkm) -

1 3 α 1 δ km , [H (2) ] km = 1 2 ( pqk lnm + pqm lnk )F (lp)nq -α 2 δ km , [H (3) ] km = 1 2 lnq ( pqm F (lp)nk + pqk F (lp)nm ), [H (4) ] km = 1 2 lnq ( pqm F (lpk)n + pqk F (lpm)n ), (11)
[H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] kmn = [S [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] (kmn) -1 5 (δ mn [a [START_REF] Nguyen | Nanoscale flexoelectricity[END_REF] ] k + δ km [a [START_REF] Nguyen | Nanoscale flexoelectricity[END_REF] ] n + δ kn [a [START_REF] Nguyen | Nanoscale flexoelectricity[END_REF] ] m ), [H [START_REF] Liu | Effect of flexoelectricity on electrostatic potential in a bent piezoelectric nanowire[END_REF] ] kmn = [S [START_REF] Liu | Effect of flexoelectricity on electrostatic potential in a bent piezoelectric nanowire[END_REF] ] (kmn) -1 5 (δ mn [a [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] k + δ km [a [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] n + δ kn [a [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] m ), (12) [H] klmn = F (klmn) -1 7 (δ kl [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] mn + δ km [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] ln + δ kn [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] lm + δ lm [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] kn + δ ln [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] km + δ mn [H [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] kl ) [S [START_REF] Bursian | Nonlocal piezoelectric effect[END_REF] ] kmn = pqn F (pk)qm , [S [START_REF] Liu | Effect of flexoelectricity on electrostatic potential in a bent piezoelectric nanowire[END_REF] ] kmn = pqn F (pkm)q .

-
(14) Next, by using the Cartan method, the second-, thirdand fourth-order harmonic tensors H (i) , H (i) and H in Eq. ( 8) can be rewritten as follows:

H (i) = α (2) 0i U 0 + α (2) 1i U 1 + β (2) 1i T 1 + α (2) 2i U 2 + β (2) 2i T 2 with i = 1, 2, 3, 4, H (i) = α (3) 0i U 0 + α (3) 1i U 1 + β (3) 1i T 1 + α (3) 2i U 2 + β (3) 2i T 2 + α (3) 3i U 3 + β (3) 3i T 3 with i = 1, 2, H = α (4) 0 U 0 + α (4) 1 U 1 + β (4) 1 T 1 + α (4) 2 U 2 + β (4) 2 T 2 + α (4) 3 U 3 + β (4) 3 T 3 + α (4) 4 U 4 + β (4) 4 T 4 , (15)
where U i , T i , U i , T i , U i and T i are the tensors involved in the Cartan decomposition, whose explicit expressions can be found in Forte and Vianello [START_REF] Forte | Symmetry classes for elasticity tensors[END_REF][START_REF] Forte | Symmetry classes and harmonic decomposition for photoelasticity tensors[END_REF] or Le Quang and He [START_REF] Quang | The number and types of all possible rotational symmetries for flexoelectric tensors[END_REF] .

With the help of the harmonic and Cartan decompositions presented above, it can be shown that the number of all possible rotational symmetry classes for all flexoelectric tensors is 12. These 12 symmetry classes are characterized by the conjugates of the following 12 sub-groups of SO(3):

I, {Z r }, {D r }, {O}, {T }, {SO(2)}, {O(2)}, SO(3) (16)
where 2 ≤ r ≤ 4.

Alternatively, the 12 sets

F(I), F(Z r ), F(D r ) with 2 ≤ r ≤ 4, F(O), F(T ), F(SO(2)), F(O(2)), F(SO(3))
form a partition of the flexoelectric tensor space F:

F = F(I) ∪ 4 r=2 F(Z r ) ∪ 4 r=2 F(D r ) ∪ F(O) ∪ F(T ) ∪ F(SO(2)) ∪ F(O(2)) ∪ F(SO(3)). ( 17
)
For more detail about the derivation of these results, the reader can refer to our previous work [START_REF] Quang | The number and types of all possible rotational symmetries for flexoelectric tensors[END_REF] .

IV. COMPACT EXPLICIT MATRIX REPRESENTATIONS OF THE FLEXOELECTRIC TENSOR

In order to obtain explicit matrix expressions of the fourth-order flexoelectric tensor for all the 12 rotational symmetry classes, we first adopt the following reduced suffix notations for the gradient E of the infinitesimal strain tensor ε and for the fourth-order flexoelectric tensor F ∈ F: where γ defined in Table II and the summation convention does not apply on j and k.

Ẽγ = √ 2(1 -δ jk ) + δ jk E jkl , (18)
Fiγ = √ 2(1 -δ jk ) + δ jk F ijkl , (19)
According to the previously reduced notation rules, the constitutive relation (3) between the electric polarization and the gradient of the infinitesimal strain tensor can now be written in the following matrix form:

p = F • Ẽ ( 20
)
where the vectors p and Ẽ are specified by

p = [p 1 p 2 p 3 ] T , Ẽ = [ Ẽ1 Ẽ2 . . . Ẽ18 ] T , ( 21
)
and F is the 3 × 18 flexoelectric matrix whose block matrix representation has the general expression

F =   D x C xy C xz J x C yx D y C yz J y C zx C zy D z J z   .
Here, the diagonal block matrices D With the previously reduced notation rules, the action of an orthogonal tensor Q ∈ O(3) on a flexoelectric tensor F can be expressed in the following simple and explicit matrix form:

Q ir Q jm Q kn Q ls F rmns = Q ir Frζ Qζγ ( 22
)
where Qζγ , standing for the components of the 18 × 18 matrix Q, are given by

Qζγ = √ 2(1 -δ jk ) + δ jk × √ 2(1 -δ mn ) + δ mn Q jm Q kn Q ls . (23)
Above, γ and ζ are defined in Table II and the summation convention does not hold for j, k, m and n. By using the procedure elaborated in Le Quang and He [START_REF] Quang | The number and types of all possible rotational symmetries for flexoelectric tensors[END_REF] to construct a flexoelectric tensor F exhibiting a required symmetry and by exploiting the harmonic and Cartan decompositions in Eqs. ( 8) and ( 15) together with the Cartan decomposition parameters α (•)

• and β III, we can exactly calculate the number of independent components contained in a flexoelectric tensor F ∈ F belonging to a given symmetry class. The corresponding results are shown in Table IV.

(•) • provided in Table
Next, by adopting the above reduced suffix notations for the fourth-order flexoelectric tensor F, we obtain the explicit matrix forms of the elementary block matrices of (D, C, J)-type as well as the explicit matrix form of F for each symmetry class. These explicit block matrix forms of F are provided as follows:

Identity class:

I FI =    D (5) x C (5)
xy C

(5)

xz J

(3) x

C

(5) yx D

(5) y

C

(5) yz J

(3) y

C

(5)

zx C

(5) zy D

(5) z

J (3) z    , (24)
Cyclic classes:

Z r FZ2 =    D (5) x C (5) xy 0 0 C (5) yx D (5) y 0 0 0 0 D (5) z J (3) z    , (25) FZ3
=    D (4) C (4) xy C
(2) xz

C

(2) xy D [START_REF] Catalan | The effect of flexoelectricity on the dielectric properties of inhomogeneously strained ferroelectric thin films[END_REF] Spatial classes:

yz • P 1 -C (4) xy D (4) C (2) yz -C (2) xz • P 1 C (1) zx C (1) zy D (3) z J (1) z    , (26)
0 0 0 0 D (3) z J (1) z    , (27)
FSO(2) =    D (4) C (4) xy 0 0 -C (4) xy D (4) 0 0 0 0 D (3) z J (1) z    , (28)
Dihedral classes: D r FD2 =    D (5) x 0 0 0 0 D (5) y 0 0 0 0 D (5) z 0    , (29) FD3
=    D (4) 0 0 C (2) yz • P 1 0 D (4) C (2) yz 0 0 C (1) zy D (3) z 0    , (30) FD4
=   D (5) 0 0 0 0 D (5) 0 0 0 0 D (3) z 0   , (31)
FO(2) =   D (4) 0 0 0 0 D (4) 0 0 0 0 D (3) z 0   , (32)
I {Z2} {D2} {Z3} {D3} {Z4} {D4} {T } a {O} b {O(2)} {SO(2)} {SO(3)} α1 • • • • • • • • • • • • α2 • • • • • • • • • • • • a (i) 1 • • • • • • • • • • • • a (i) 2 • • • • • • • • • • • • a (i) 3 • • • • • • • • • • • • α (2) 0i • • • • • • • • • • • • α (2) 1i • • • • • • • • • • • • β (2) 1i • • • • • • • • • • • • α (2) 2i • • • • • • • • • • • • β (2) 2i • • • • • • • • • • • • α (3) 0i • • • • • • • • • • • • α (3) 1i • • • • • • • • • • • • β (3) 1i • • • • • • • • • • • • α (3) 2i • • • • • • • • • • • • β (3) 2i • • • • • • • • • • • • α (3) 3i • • • • • • • • • • • • β (3) 3i • • • • • • • • • • • • α (4) 0 • • • • • • • • • • • • α (4) 1 • • • • • • • • • • • • β (4) 1 • • • • • • • • • • • • α (4) 2 • • • • • • • • • • • • β (4) 2 • • • • • • • • • • • • α (4) 3 • • • • • • • • • • • • β (4) 3 • • • • • • • • • • • • α (4) 4 • • • • • • • • • • • • β (4) 4 • • • • • • • • • • • • a α (4) 4 = 5α
FT =   D (5) 0 0 0 0 D (5) • P 2 0 0 0 0 D (5) 0   , (33)
FO =   D (3) 0 0 0 0 D (3) 0 0 0 0 D (3) 0   , (34)
FSO(3) =   D (2) 0 0 0 0 D (2) 0 0 0 0 D (2) 0   . (35)
In Eqs. ( 25)-( 35), the superscript of each elementary block matrix of D, C or J-type denotes the number of independent material parameters; P 1 and P 2 are two cor-rection matrices. By introducing two material parameters defined as

η = 1 √ 2 ( F11 -F12 ), θ = 1 √ 2 ( F22 -F21 ), (36)
the elementary block matrices of (D, C, J)-types as well as of the correction matrices P 1 and P 2 are explicitly expressed as follows: D-type elements:

D (5) = d 1 d 2 d 3 d 4 d 5 , (37)
D (4) = d 1 d 2 η d 4 d 5 , (38)
D (3) = d 1 d 2 d 3 d 2 d 3 , (39)
D (2) = d 1 d 2 η d 2 η , (40)
C-type elements:

C (5) = c 1 c 2 c 3 c 4 c 5 , (41)
C (4) = c 1 c 2 θ c 4 c 5 , (42)
C (2) = 0 c 2 c 3 -c 2 -c 3 , (43)
C (1) = c 1 -c 1 -c 1 0 0 , (44)
J-type elements:

J (3) = j 1 j 2 j 3 , J (1) = 0 j 2 -j 2 , (45)
Correction elements:

P 1 =      0 0 0 1 0 0 0 1 1 0 0 0 0 0 1      , P 2 =      1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0      . ( 46
)
Finally, the compact explicit matrix expressions for the 12 symmetry classes of flexoelectric tensors are shown in Table V.

V. REFLECTION SYMMETRY CLASSES OF THE FLEXOELECTRIC TENSOR

First, let us introduce P(n) ∈ O(3)\SO( 3), the reflection through the plane P(n) = {x ∈R 3 | x • n = 0} perpendicular to a unit vector n, by

P(n) = I -2n ⊗ n. (47)
It can be seen from ( 47) that P is an even function of n in the sense that P(-n) = P(n). Thus, two unit normal vectors are associated to each reflection. Next, we denote by P F the set of reflection symmetry elements of F ∈ F. The reflection symmetry class of F ∈ F, symbolized by {P F }, is defined as the collection of all the conjugates of P F , namely

{P F } = {P F = RP F R T | R ∈ SO(3)}. ( 48
)
Note that {P F } represents the reflection symmetry class and not the symmetry group of F ∈ F. We introduce the unit vector

r i (θ) = sin θe j + cos θe k (49)
with {i, j, k} being a cycle permutation of {1, 2, 3}, and the following sets of reflection transformations:

• P h being the set which contains only the reflection P(e 3 ). The label h means "horizontal ".

• P v k being the set defined by P v k = {P(r 3 ( 2pπ k )} 1≤p≤k and consisting of k elements. The label v means "vertical ".

• P hv k (k ≥ 1) being the set defined by the k elements of P v k completed by P(e 3 ): this set comprises k +1 elements.

• P O designating the cubic set of nine reflections with respect to the nine planes of which the normals of 6 pass through the center of each edge of a regular cube and the normals of 3 through the center of each face of the latter.

• P I denoting the icosahedral set consisting of fifteen reflections with respect to the fifteen planes whose normals pass through the center of each edge of a regular icosahedron;

• P O(3) representing the set composed of all reflections P(n) with n ∈ S 2 where S 2 is the unit sphere defined by

S 2 = {x ∈ R 3 | x = 1}.
We can show that the space F of flexoelectric tensors is divided into the 8 reflection symmetry classes which are characterized by the following 8 sets of reflection transformations:

∅, {P h }, {P hv2 }, {P v3 }, {P hv4 }, {P hv∞ }, {P O }, {P O(3) }.
(50) The characteristics of each reflection symmetry class, and especially its link with the rotational ones, are detailed in Table VI. It is important and interesting to remark that: (i) unlike the results as obtained by Chadwick et al. [START_REF] Chadwick | A new proof that the number of linear anisotropic elastic symmetries is eight[END_REF] for the space of fourth-order elasticity tensors, according to which the classifications by rotational symmetry groups and by reflection symmetry planes give the same response, the number of rotational symmetry classes for the space of fourth-order flexoelectric tensors is 12 while the one of reflection symmetry classes is only 8; (ii) these 8 reflection symmetry classes for the space of flexoelectric tensors are exactly identical to the ones for the space of elasticity tensors, even through a fourth-order flexoelectric tensor is algebraically more complex than a fourthorder elasticity tensor.

VI. A GRAPHIC METHOD FOR IDENTIFYING REFLECTION AND ROTATIONAL SYMMETRIES

The important question now arises as to how to identify the reflection symmetry and rotational symmetry and of a given flexoelectric tensor from the knowledge of its matrix relative to a basis. The present section aims at elaborating a simple but efficient graphic method to answer this question.

Francois et al. [START_REF] Francois | Determination of the symmetries of an experimentally determined stiffness tensor: Application to acoustic measurements[END_REF] initiated a graphic approach to identifying the reflection symmetry planes that a given fourthorder elastic tensor has. This approach is based on the notion of "pole figures". Owing to the fact that the reflection symmetry classes of the fourth-order elastic tensor are identical to its rotational symmetry classes, the identification of the formers leads also the one of the latters. However, when the fourth-order elastic tensor is concerned, the situation is much more complicated, since its reflection symmetry classes are 8 while its reflection symmetry classes are 12. Thus, in this section we extend TABLE V. Number of independent material parameters and compact explicit matrix expression for each of the 12 symmetry classes of flexoelectric tensors classes of flexoelectric tensors

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !"#$%#&#$'()*+ ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! "#$#%&'$'%()* +*,- ! ! ! ! ! ! ! ! ! ! ! ! !"#$%&'()*+ ,-./ ! ! ! % ! ! ! ! ! 0 % % ! ! % % % ! ! 0 ! ! ! ! ! % " % ! " ! ! % " ! ! ! ! " " " 1 " % % % % ! % ! % % ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! "#$%&$#&'() *+ ,-./ ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! "#$%&'()*+, -./0 ! ! ! ! 1 ! ! ! ! 1 ! ! ! ! & ! ! ! ! ! ! ! ! ! ! ! ! & ! & ! & ! ! ! ! "#$%&'!(&) *+ ,-+. ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !"#$% ! &'( ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !"#$%&'(%) *+ ,-. ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! "#$%&'#(%&) ! *+, ! ! ! ! ! ! ! ! x x ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ Trans. Iso. O(2) (7) • • • • • • • • ⋅ ⋅ ⋅ • • • ⋅ ⋅ ⋅ ⋅ • • ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ !"#$%&'(%)&'*+,-.',/-. ! ! ! 0 ) ) ) ! 0 ! ! ! ! " " " ! ! ! " " " 1 "
) ! ! " " " " " " " " " " " " " " " " " " " " the graphic approach of Francois et al. [START_REF] Francois | Determination of the symmetries of an experimentally determined stiffness tensor: Application to acoustic measurements[END_REF] to being able to identify not only the reflection symmetry of a given flexoelectric tensor but also its rotational symmetry.

x x x ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ Isotropic SO(3) (2) • • • • x • • • ⋅ ⋅ ⋅ • x • ⋅ ⋅ ⋅ ⋅ x ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
Let n ∈ S 2 be the unit vector relative to the reflection transformation P(n) = I -2n ⊗ n through the plane P(n). With no loss of generality, n is expressed by n = sin θ cos φe 1 + sin θ sin φe 2 + cos θe 3 (51)

where (φ, θ) ∈ [0, 2π[×[0, π[ denote, respectively, the longitude and colatitude angles relative to a system of spherical coordinates. Then, by considering, for a given flexoelectric tensor F ∈ F, the function:

L(θ, φ) = P(θ, φ) * F -F = P(θ, φ) • F • P(θ, φ) -F (52
) in which • is the Frobenius norm inherited from the scalar product on F, P(θ, φ) is a reflection operator parametrized with the longitude and colatitude angles and P(θ, φ) is a 18 × 18 matrix whose components are obtained by replacing Q with P in (23). Finally, the vanishing loci of L(θ, φ) give the unit normals to the symmetry planes that F has.

Concretely, the function ( 52) is numerically evaluated in a discrete way. Precisely, we introduce

M ij = L(θ i , φ j ) with θ i = i 2π N and φ j = j 2π N
where the number N depends on the degree of numerical accuracy required. In our computations, N is set to be equal to 160. To evaluate the function L(θ, φ), we first use the matrix representations of flexoelectric tensors presented in Section IV. In addition, the numerical values of the components of a flexoelectric tensor are determined as random integers picked-up in the range {-10, 10}. We show, in Figure 1, the loci of the zeros of L(θ, φ) plotted on the θ -φ plane for all rotation symmetry classes and all reflection symmetry classes. It can be seen from Figure 1 that the number of symmetry planes for a given rotation symmetry class or reflection symmetry class coincides exactly with the one provided in Table VI. This also constitutes a validity verification of our theoretical results.

In addition, it can be observed from Figure 1 that, even through the flexoelectric tensors belonging to both rotational symmetry classes {Z 3 } and I do not exhibit any reflection symmetry plane, we can differentiate them since the loci of the function L(θ, φ) for the flexoelectric tensors belonging to {Z 3 } are periodic in the φ-direction with period 2π/3 while the ones for the flexoelectric tensors appertaining to I are not periodic. Similarly, even if the flexoelectric tensors belonging to the rotational symmetry classes {D 2 } and T possess the same number of reflection symmetry planes, the loci of the function L(θ, φ) for flexoelectric tensors belonging to {T } are periodic in the φ-direction with period π/2 but the counterpart of the flexoelectric tensors appertaining to {D 2 } are periodic in φ-direction with period π.

The foregoing graphical approach is now applied to a given flexoelectric tensor F ∈ F with a given angle ψ ∈ [0; 2π[ for identifying all invariant directions defined by n = sin θ cos φe 1 + sin θ sin φe 2 + cos θe 3 in the sense that F is unchangeable under the rotational transformation action Q(n, ψ). As before, by using the well-known Rodrigues expression of Q(n, ψ), i.e.

Q(n, ψ) = cos(ψ)I -sin(ψ) • n + [1 -cos(ψ)]n ⊗ n (53)
in which denotes the Levi-Civita third-order tensor, and by introducing the following function

J(θ, φ, ψ) = Q(θ, φ, ψ) * F -F = Q(θ, φ, ψ) • F • Q(θ, φ, ψ) -F (54)
where Q(θ, φ, ψ) is a rotation operator parametrized with the longitude, colatitude and rotation angles; Q(θ, φ, ψ) is a 18 × 18 matrix whose components are defined by (23), the vanishing loci of J(θ, φ, ψ) allow us to obtain the invariant directions n that F possesses.

We illustrate, in Figure 2, the vanishing loci of the function J(θ, φ, 2π/3) plotted on the θ -φ plane for a flexoelectric tensor F belonging to the rotational symmetry class {Z 3 } whose matrix representation is provided in Section IV. It can be seen from Figure 2 that there exists only one invariant axis spanned by n = e 3 or n = -e 3 . This is in agreement with the fact that the rotational symmetry class {Z 3 } contains Q(e 3 , 2π/3).

By combining the knowledge of the matrix representations of all rotational symmetry classes with pole figures, we can finally identify the rotational symmetry class to which a given flexoelectric tensor belongs. The corresponding identification procedure is summarized in Figure 3.

VII. CONCLUDING REMARKS

Flexoelectricity is an electromechanical phenomenon which has a great number of potential applications including energy harvesting, sensors, actuators and biotechnology. A full understanding of the fourth-order flexoelectric tensor is essential not only to the fundamental theory of flexoelectricity but also to all possible applications of flexoelectricity. In the present work, which may be viewed as a continuation of our previous one [START_REF] Quang | The number and types of all possible rotational symmetries for flexoelectric tensors[END_REF] , compact explicit matrix representations of the flexoelectric tensor have been provided for all the 12 possible rotational symmetry classes, so as to facilitate its use in various situations; the reflection symmetry classes of the flexoelectric tensor have been also determined and shown to be identical to those of the fourth-order elastic tensor; a simple and efficient graphic method for identifying the rotational symmetry and reflection symmetry of a given flexoelectric tensor has been elaborated and illustrated. These results contribute to developing the continuum theory of flexoelectricity and rendering the use of this theory easier in various anisotropic cases. The matrix FABO3 has 3 independent components, i.e., the longitudinal flexoelectric coefficient F 1111 , transverse flexoelectric coefficient F 2112 and shear flexoelectric coefficient F 2121 . The values of F 1111 , F 2112 and F 2121 determined experimentally and computationally can be found in Wang et al. [START_REF] Wang | Flexoelectricity in solids: Progress, challenges, and perspectives[END_REF] and Shu et al. [START_REF] Shu | Flexoelectric materials and their related applications: A focused review[END_REF] for some Pervoskites.

In parallel with the formulation used in the present work in which the electric polarization vector p is linearly related to the strain gradient E through the fourth-order flexoelectric tensor F, there is another formulation of the electric polarization vector p linearly connected to the second-order derivative of the displacement vector, U = ∇∇u, as follows (see e.g. Hong et al. [START_REF] Hong | First-principles theory and calculation of flexoelectricity[END_REF] ):

p i = F I ijkl U jkl or p i = F I ijkl u j,kl . (55)
Here F I ijkl , the tensor component of the type-I flexoelectric tensor F I , possesses the index permutation symmetry F I ijkl = F I ijlk . Compared with F, called also type-II flexoelectric tensor, more suitable not only for formulating the thermodynamic theory of ferroelectric materials but also for making comparisons with experimental measurements, the definition of the type-I flexoelectric tensor F I is complicated for mathematical derivations of the microscopic theory of flexoelectricity. It can be shown that the connections between F I and F are given by

F I ijkl = 1 2 (F ijkl + F ijlk ), F ijkl = F I ijkl + F I ikjl -F I iljk .
(56) In addition, due to the fact that both flexoelectric tensors F I and F exhibit mathematically the same permutation symmetry with respect to two indexes, the flexoelectric tensors F I and F will possess the same rotation and reflection symmetry classifications. Moreover, by adopting can be expressed in a 3 × 18 matrix form FI with components F I iγ . Consequently, the number of independent material parameters and compact explicit matrix expression for each of the 12 symmetry classes of flexoelectric tensors F I are exactly identical to the ones of F.

Finally, it is interesting and important to remark that the methods and results presented in the present work are directly applicable to the fourth-order flexomagnetic tensor [START_REF] Lukashev | Flexomagnetic effect in frustrated triangular magnetic structures[END_REF] .

  x , D y and D z of Dtype and extra-diagonal block matrices C xy , C xz , C yx , C yz C zx and C zy of C-type have the same size 1 × 5 while the extra-diagonal block matrices J x , J y and J z of J-type are of the size 1 × 3.

FIG. 1 .FFF

 1 FIG. 1. Loci of the zeros of the function L(θ, φ) plotted on the θ -φ plane for all rotation symmetry classes and all reflection symmetry classes.

FIG. 3 .

 3 FIG. 3. Procedure for identifying the rotational symmetry class to which a flexoelectric tensor belongs with respect to a given basis

TABLE I

 I

		. List of notations
	Notations	Descriptions
	p, pi	Electric polarization vector and its components
	F, F ijkl	Fourth-order (type-II) flexoelectric tensor and its components
	F, Fiα	Matrix representation of the flexoelectric tensor and its components
	F I , F I ijkl	Fourth-order type-I flexoelectric tensor and its components
	D, D ijk	Third-order piezoelectric tensor and its components
	ε, εij	Infinitesimal strain tensor and its components
	E, E ijk	Strain gradient tensor and its components
	Ẽ, Ẽα	Vector representation of the strain-gradient tensor and its components
	H, H, H	Second-, third-, fourth-order harmonic tensors
	F	Space of flexoelectricity tensors
	Q	Generic orthogonal tensor
	Q(a, θ)	Rotation about a through an angle θ
	G(F)	Rotation symmetry group of F
	{Gi}	Rotation symmetry class
	Zn (n ≥ 2)	Cyclic group of order n, generated by the n-fold rotation Q (e3, θ = 2π/n)
	Dn (n ≥ 2)	Dihedral group of order 2n generated by Zn and Q(e1, π)
	SO(3)	Tridimensional rotational group
	O(3)	Tridimensional orthogonal group
	T	Tetrahedral group of 12 elements generated by Q

  α 1 15 (δ kl δ mn + δ km δ ln + δ kn δ lm ). (13) Above, δ ij and ijk are the Kronecker delta and permutation symbol, respectively; either • i1i2...in or [•] i1i2...in is the component of an nth-order tensor •; • (i1i2...ir)ir+1...in denotes the average of r! components obtained by permuting the indices i 1 , i 2 , . . . , i r in all possible ways; the third-order tensors S 1 and S 2 are defined as

TABLE II .

 II Suffix notation correspondences between (j, k, l) and γ.

	(j, k, l) or (m, n, s)	γ or ζ
	(1,1,1)	1
	(2,2,1)	2
	(1,2,2) or (2,1,2)	3
	(3,3,1)	4
	(1,3,3) or (3,1,3)	5
	(2,2,2)	6
	(1,1,2)	7
	(1,2,1) or (2,1,1)	8
	(3,3,2)	9
	(2,3,3) or (3,2,3)	10
	(3,3,3)	11
	(1,1,3)	12
	(1,3,1) or (3,1,1)	13
	(2,2,3)	14
	(2,3,2) or (3,2,2)	15
	(1,2,3) or (2,1,3)	16
	(1,3,2) or (3,1,2)	17
	(2,3,1) or (3,2,1)	18

TABLE III .

 III Zero components • and non-zero independent material parameters • contained in a flexoelectric tensor belonging to a given symmetry class Parameters

TABLE IV .

 IV Number of independent material parameters contained in a flexoelectric tensor belonging to a given symmetry class

	Symmetry class	I {Z2} {D2} {Z3} {D3} {Z4} {D4} {T } {O} {O(2)} {SO(2)} {SO(3)}
	Number of independent components 54 28	15	18	10	14	8	5	3	7	12	2

TABLE VI

 VI

	. Symmetry plane stratification of the space F of fourth-order flexoelectric tensors
	System	Reflection	Number of reflection	Rotational
		symmetry class symmetry planes	symmetry class
	Triclinic	∅.	0	I, {Z3}
	Monoclinic	{P h }	1	{Z2}, {Z4}, {SO(2)}
	Orthotropic	{P hv 2 }	3	{D2}, {T }
	Trigonal	{Pv 3 }	3	{D3}
	Tetragonal	{P hv 4 }	5	{D4}
	Trans. isotropic	{P hv∞ }	∞ + 1	{O(2)}
	Cubic	{PO}	9	{O}
	Isotropic	{P O(3) }	∞ 3	{SO(3)}

TABLE VII .

 VII Suffix notation correspondences between (j, k, l) and γ. the reduced suffix notations described in Table VII the fourth-order flexoelectric tenor F I with component F I ijkl

	(j, k, l)	γ
	(1,1,1)	1
	(1,2,2)	2
	(2,1,2) or (2,,2,1)	3
	(1,3,3)	4
	(3,1,3) or (3,3,1)	5
	(2,2,2)	6
	(2,1,1)	7
	(1,1,2) or (1,2,1)	8
	(2,3,3)	9
	(3,2,3) or (3,3,2)	10
	(3,3,3)	11
	(3,1,1)	12
	(1,1,3) or (1,3,1)	13
	(3,2,2)	14
	(2,2,3) or (2,3,2)	15
	(3,1,2) or (3,2,1)	16
	(2,1,3) or (2,3,1)	17
	(1,2,3) or (1,3,2)	18

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.