Baptiste Ristagno
email: baptiste.ristagno@univ-lorraine.fr

Dominique Giraud

Julien Fontchastagner

Denis Netter

Noureddine Takorabet

Geoffrey Devornique

Nicolas Labbe

FEM using projection of physical properties suitable for movement modeling and optimization processes

Keywords: Finite Element Method, Mesh, Projection 1

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

FEM using projection of physical properties suitable for movement modeling and optimization processes.

1 Introduction

Nowadays, Finite Element Method (FEM) is widely used in electromagnetic field modeling which provides very accurate results compared to experimentations. It also allows the description of complex geometries and problems with a very simple implementation. These advantages explain the spread of FEM in many fields. Nevertheless, for complex three dimensional problems, it requires huge computational resources and consequently very long CPU time. Optimization processes usually need a high number of calls of the objective function. Therefore, FEM does not seem suitable in 3D optimization processes.

In literature, several modeling methods, such as hybrid models which combine finite element and reluctance network, have been suggested in order to reduce simulation time [START_REF] Devornique | Hybrid model: Permance network and 3-d finite element for modeling claw-pole synchronous machines[END_REF]. In this article, the authors propose a method allowing meshing step saving in any iterative process such as movement modeling. The mesh being fixed, all physical parameters (properties, supply or armature movement) are implemented by projecting mathematical functions on the mesh grid. This method is applied to a usual permanent magnet synchronous machine (PMSM) in 2D and in 3D. Interesting results are obtained and discussed.

Methodology

2.1 A brief overview of the usual way of using FEM

In the FEM method, each elementary part of the initial geometry is implemented in the software. This discretization is then used to build an adapted mesh in which all nodes match with material boundaries. Behavior law and physical properties of all materials are defined based on this mesh. Consequently, we can take the example of the reluctivity ν.

It is a piecewise function defined as ν i with i corresponding to different material regions Ω i . The same characteristic can be highlighted on the current density j or the remanent flux density b r , both often used in magnetic problems, as seen in Figure 1.

The global weak formulation is applied to the whole domain using all the piecewise functions as explained by [START_REF] Bastos | Magnetic Materials and 3D Finite Element Modeling[END_REF]. For example, a permanent magnet synchronous machine implies the weak formulation along with vector potential a seen in (1). where Ω is the whole domain (∪ Ω i , i ∈ [1; number of subdomains], Ω l magnets subdomains with l ∈ [1; number of magnet subdomains] and Ω m source subdomains with m ∈ [1; number of source subdomains]. Numerical integration method applied to weak formulation is strongly linked to spatial discretization of geometry. Some iterative processes as optimization problem or movement modeling need to modify the geometry and consequently the weak formulation support. As iterative processes required a huge number of steps, all these pre-processing part reevaluations represent a non negligible waste of time, especially in 3D.

k Ω k ν k ( curl a ) curl a • curl a dΩ k - l Ω l ν l ( curl a ) b rl • curl a dΩ l - m Ωm j m • a dΩ m = 0, ∀a ∈ H 0 (curl, Ω) (1) Physical properties definition

Proposed approach

The authors propose to develop a simple method to avoid the pre-processing part reevaluations in an iterative process, similar to the moving band method in the 2D rotating movement modeling. In the case of optimization processes only few methods exist to reduce the time of each evaluation as in [START_REF] Devornique | Three dimensional pole shape optimization of claw pole machines based on a hybrid model[END_REF] in which efficient FE models are developed to find a balance between precision and CPU time. That is why the authors proposed, by dint of a trick, a different approach in the use of FEM which mathematically respects the standard weak form. Two crucial points have to be pointed out concerning the traditional way of using FEM :

• The discretization of material boundaries is a good basis for the mesh. Indeed, all nodes and edges perfectly match with boundaries, allowing a good evaluation of behavior law and physical properties at each node, thus at each boundaries. In conclusion, discretization is useful to correctly define material boundaries.

• According to the previous spatial discretization, material properties are implemented in each spatial area, defining a piecewise function of the global domain.

The implementation of piecewise functions (an example is given in ( 2)) has to be respected. On the other hand, the fact that nodes do not match with boundaries only implies some imprecisions. Considering these facts, the authors chose to define material properties by projection of mathematical functions on a mesh grid to represent the geometry. In other words, the first step is to only define the external boundary of the problem. Using these boundaries, a fixed mesh can be created and used as a support for physical property functions. Indeed, physical properties are no more implemented in discrete regions, as seen in Figure 1, but projected on the whole meshed domain by using space dependent piecewise functions as seen in Figure 2. The whole geometry and supply can be defined thanks to this method. br(x(t),y(t),z(t)) j(x(t),y(t),z(t)) ν(║b ,x(t),y(t),z(t)) ║

Figure 2: Principle of the proposed method applied on a conventional PMSM

In the case of the previous PMSM, presented in Figure 1 and Figure 2, two physical property functions, the reluctivity ν(x, y, z, t, b ) and the remanent induction b r (x, y, z, t), along with a source, the current density j(x, y, z, t), are required to describe the complete problem. In the example of the current density, it can be described as in (2), where R 1 and R 2 are inner and outer radii of slots, Θ the angles origin, τ r the slot pitch, ω the pulsation of slot current, I max the maximum slot current and S slot the slot surface.

j(x, y, z, t) =                                      if x 2 + y 2 ≥ R 1 & x 2 + y 2 ≤ R 2 for k = 1 : 1 : Number of slots if atan y x > Θ + 2kτ r & (atan y x < Θ + τ r + 2kτ r Imax S slot cos(ωt + k 2π 3 ) e z else 0 e z endif endfor else 0 e z endif (2)
With this approach, defining the external boundary and an adapted fine mesh (discussed in 2.3) is enough to build the model. The mesh is kept unchanged for all iterations of any iterative process, which saves CPU time by avoiding pre-processing evaluatation step. Geometry is defined by projection of physical property functions mentioned above, as seen in Figure 2. Mathematically, the weak formulation seen in ( 3) is kept unchanged as in (1). Obviously, both formulations are the same, there are only seen from a different point of view.

ν(x(t), y(t), z(t), curl a ) curl a • curl a dΩ - Ω ν(x(t), y(t), z(t), curl a ) b r (x(t), y(t), z(t)) • curl a dΩ - Ω j(x(t), y(t), z(t)) • a dΩ = 0, ∀a ∈ H 0 (curl, Ω) (3)
Basically, the aforementioned method has some similarities with the mesh-less method presented by [START_REF] Belytschko | Element-free galerkin methods[END_REF] although finite elements are still used to lead to well conditioned matrices. Furthermore, the geometry of finite elements in the proposed method is standard and completely decoupled from the geometry (computer aided design (CAD)) of the simulated device. Whereas in the meshless method, the geometry of the elements is arbitrary (no elements, only nodes) but node placements are strongly linked to the geometry of the device.

Drawbacks and discussions

In this section, the principal drawbacks of the method will be discussed and some advantages will be highlighted. The major drawback is directly linked to the spatial discretization of the usual way of using FEM. As previously mentioned, the interest of the spatial geometry discretization is every node of mesh corresponds to material boundaries. The fact that the material properties are evaluated each node of mesh allows for very well-defined material boundaries. For instance, a crucial boundary for rotating machines is the air-gap which small variation of thickness deeply impacts both global and local variables.

To avoid re-meshing step, the spatial discretization is replaced by the projection of physical properties on a mesh grid. In this condition, mesh nodes do not match anymore with material boundaries. It gets harder to implement clear and perfect boundaries. Indeed, between two nodes, the mathematical piecewise function is defined by inducing an evaluatation of physical properties thanks to a linear interpolation along the edge [START_REF] Bossavit | Edge-element computation of the force field in deformable bodies[END_REF]. This situation leads to irrregular boundaries as seen in Figure 3. The boundary is defined between two nodes, so that it can lead to a deviation of thickness of one element. It is interesting to point out the fact that a ruled mesh provides more regular boundaries as seen in Figure 4, but there is still a deviation of thickness. A first solution to reduce the deviation is to reduce the size of the elements to converge to the real boundary. That is to say, an adapted fine mesh can be used to limit the impact of the linear interpolation deviation. Nevertheless, adapting the mesh grid to the smaller mesh element size is a nonsense and will raise the CPU time compared with the usual method because of the high number of Dofs.

Thus, the authors propose to adapt the mesh grid defining some spatial discretizations whose only purpose is to allow a control on the mesh. These geometrical elements are not correlated with real geometry but can be placed wherever the users want to set a fine mesh. Indeed, in the pre-processing part in which external boundaries of the problem are defined (so before the projection step), some non physical points or lines can be placed in interesting areas. Those geometrical elements allow for implementing mesh size near points or number of nodes on lines. Taking the example of Figure 3, the size of the mesh element of the whole domain is about 1 p.u., leading to a deviation of the boundary of 1 p.u. compared to the expected boundary. Defining a line near the boundary allows for reducing the mesh element size only close to the boundary while letting it wide elsewhere. For instance, the mesh element size is reduced by a factor 3 on Figure 5, leading to a more precise evaluation of the boundary, with a deviation of 0.33 p.u. only. Obviously, it is more time-consuming and the meshing time between Figure 3 and Figure 5 is increased by a factor 2. This technique can be illustrated on an electrical machine: the boundary at the air-gap of a PMSM is clearly important, which is why the mesh has to be controlled, but a low deviation at the external boundary of the yoke will not have a deep impact on performance. According to the previous paragraph, there is no mathematical compensation to the nodal linear interpolation. Resulting irregular boundaries can be seen as rough surfaces. This lies in the fact that the used software evaluates pre-processing properties and sources only on the mesh nodes. The linear interpolation between two nodes is quite unavoidable in this case, nevertheless it is imaginable to evaluate properties and sources along mesh edges with a higher order of interpolation of pre-processing functions using maybe another software.

In addition to the pre-processing time saving, the proposed method is mathematically interesting. The iterative process of resolution (the Newton-Raphson method for example) is better initialized than in a multi-static resolution. Indeed, since the mesh is kept unchanged, the solution of non-linear timestep n can be easily used as initialization of non-linear timestep n + 1. So the CPU time saving is expected to be higher than the saving of pre-processing part.

Tests on an application

In order to compare both methods, the same software is used : Onelab interface based on the finite element mesh generator Gmsh [START_REF] Geuzaine | Gmsh: A three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF] and the general finite element solver GetDP [START_REF] Dular | Getdp reference manual: the documentation for getdp, a general environment for the treatment of discrete problems[END_REF]. It was chosen for its flexibility with weak formulations and physical property definitions. All simulations are run with vector magnetic unknowns and magnetostatic formulation.

Permanent Magnet Synchronous Machine

The study of a usual permanent magnet synchronous machine, 12 slots and 4 poles, is carried out. The first method lies in drawing geometry by CAD and then implement all physical properties in every subdomain created, as seen in Figure 1. The same machine is then implemented by projection of physical functions, as seen in Figure 2. In both cases, particular attention is paid to implement a fine mesh near the airgap. Indeed, the rotation step and mesh element size are equalized and three nodes are placed in radial length of the airgap. In the first method, 15 nodes are implemented on a half slot pitch whereas in the second one, 360 nodes are implemented on the circumference. That is to say, both models have 360 nodes on their circumference with a one-degree rotation step. But the second model implies a constant mesh element size, as in the first one, the number of nodes on a half slot pitch is constant.

Flux density maps of both machines, in magnetostatics 2D and 3D, are compared first. Qualitative correspondence are observed on Figure 6 and Figure 7, which represents an encouraging first step for the validation of the proposed method. Nevertheless, an iterative process has to be carried out to check the expected advantages of projection method and a complete statistical study is required to try to avoid any confirmation bias. Comparisons are made between the proposed method and a basic approach without any special technique (as a moving band or a sliding surface) which can be considered as a "naive" method. In the case of basic movement modeling, a pre-processing part (geometry and mesh) is performed at each position and the system is solved. For a number of positions N b positions , the system is meshed and solved N b positions times, so:

CP U T ime = N b positions * (M eshing T ime + Resolution T ime )
In the case of the proposed movement modeling, the meshing step is avoided at each iteration. Moreover, due to the low complexity of the fixed mesh grid, the meshing step is faster than the usual method for same Dofs number. This saved time is hard to evaluate because it depends of the problem geometry but it can be represented by the ε variable (∈]0; 1]) in the following. The non-linear convergence algorithm being also better initialized in the case of the proposed method, the resolution time is reduced by a factor depending on geometry and supply. CPU time with expected time saving are sum up in Table 1.

Comparison of 3D models

Proposed method Classical method Multistatic torque is computed for both models in 2D and 3D as it can be seen on Figure 8 and Figure 9. The relative deviation between black and red curves is kept under 1% for 2D case. But as expected, the observed saved time is very low because of the unsignificant meshing time. 2D results being encouraging, 3D models are developed in order to have same number of Dofs (around 130,000). The CPU time is divided by a factor 1.7 for our proposed method, keeping a relative deviation under 5% using a Intel Core i7 2.70GHz with 32Gb of RAM.

To implement the proposed model in an optimization process, a statistical validation is needed. 500 randomized 3D machines are generated. 4 geometrical parameters vary in a limited range given in Table 2. Airgap, external radius and bore radius remain constants. For each machine, the average torque of both methods is computed and compared. Results are shown in Figure 10 for 130,000 Dofs first.

The average relative deviation is around 13% with a standard deviation of 17%. It is interesting to highlight the fact that deviations are centered on a positive value, which Figure 10: Robustness of proposed method for 130,000 Dofs means that the proposed method tends to underestimate the "usual torque". Considering average and standard deviation, 75% of the simulated machines have a relative deviation below 25% in absolute value. Moreover, it can be observed that refining mesh tends to change the histogram as it can be seen in Figure 11.

Knowing the mean value, the standard deviation and the population size, the confidence interval at 95% on mean value can be estimated. This interval is represented on Figure 12. As expected, refining mesh seems to reduce the relative deviation between both methods by improving the boundary implementation.

Concerning the ε value mentioned in Table 1, the same statistical analysis can be carried out on the 500 machines sample both for 130,000 Dofs and 260,000 Dofs. As a reminder, ε is defined in (4).

ε = M eshing T ime,projection M eshing T ime,usual (4)
It turns out that the average is 0.132 with a standard deviation of 0.003 for 130,000 Dofs and the average is 0.239 with a standard deviation of 0.005 for 260,000 Dofs. That is to say, meshing time can almost be 10 times faster with projection method than with traditional method but it seems to decrease with Dofs number. However, saved time represented by ε is strongly linked to the geometry discretization (CAD). The analysis of ε only make sense for this particular geometry with our particular choice of discretization. Nevertheless, according to previous results, the proposed model seems to be suitable for 3D optimization processes. According to [START_REF] Räisänen | How can we deal with moving objects on a fixed mesh?[END_REF], there are some additional advantages using a fixed mesh in optimization processes, such as: • Mesh generation might be difficult or even fail for certain configurations;

• Changing the mesh influences the discretization error. For instance, if optimization algorithms are applied, the target function will be polluted with re-meshing noise.

Optimization problem

Previous results show that multistatic movement of rotating machine is well-modeled by the proposed method. The authors wanted to study the efficiency of the proposed method in some optimization processes. As 2D projected models do not present a real interest, a 3D PMSM geometry optimization will be carried out. The same models as robustness study are used in this section, that is to say the 130,000 Dofs models for time saving purpose. Geometry optimization is applied to the vector solution presented in Table 2.

The process consists of a torque ripple ∆ Γ minimisation with minimum mean torque reference Γ mean,ref and maximum current density reference J max,ref .

min

X∈R (∆ Γ ) under c 1 = Γ mean,ref -Γ mean (X) ≤ 0 c 2 = J max -J max,ref ≤ 0
Torque ripple ∆ Γ is defined according to (5) in which Γ max is the maximum of instaneaous torque, Γ min is the minimum of instaneaous torque and Γ mean is the mean value of instaneaous torque.

∆ Γ = Γ max -Γ min Γ mean (5)
This optimization problem is applied to the usual model and to the proposed model as seen in the previous section. NOMAD software [START_REF] Audet | Nomad user guide[END_REF] is used to run these optimizations linked to the Onelab interface to solve the magnetic problem. NOMAD is based on a mesh adaptive direct search (MADS) algorithm and is useful for a few variables and a time consuming blackbox evaluation. Both optimization models have been identically parametrized with the same algorithm, objective, constraints and solution vector. Models are investigated by means of three problems, presented in Table 3, in which convergence and stop criteria are modified. Indeed, Problem I is a usual case of optimization process which stops when the relative deviation on the objective function is lower than a fixed limit Deviation min . Problems II and III represent two optimization processes in a given maximum time and in a given maximum number of blackbox evaluations.

Results of optimization processes are presented in Table 4. It is very interesting to point out some observations. Problem I shows that both methods converge to the same solution. Compared to the usual method, time saving is observed and the amount of blackbox evaluations increased in the proposed method. Conclusion is that one blackbox evaluation is made twice as fast as the usual method. Thus, the domain is more widely explored in the same period of time with the projection method. Problem II and III

∆ Γ ∆ Γ Constraints c 1 , c 2 c 1 , c 2 c 1 , c 2 Convergence Deviation M in - - Stop criterion - Time M ax BB Evaluation M ax
confirm this fact, the proposed method can be a little more efficient in a time limited or a blackbox evaluation limited optimization. Both methods converge to a different vector solution and results are re-calculated with the complete model.

Coupled-circuit simulation

Optimizations on magnetostatic resolution show the efficiency of the proposed method. Nevertheless, transient-state simulation are often useful for machine design, which is why the authors wanted to include in the proposed model the possibility of taking backelectromotive force (EMF) into account. Usually, the magnetic flux is computed by integration on the correct surface, from spatial discretization, using vector magnetic potential a:

ϕ i = W ire b • dS = W ire rot a • dS = slot i a • dl = 1 S slot i ( a • dl)dS
, with i ∈ [1, .., N phase ] (6)

As the proposed model has no surface defined, the normal vector to the wire crosssection dl (tangential to current density) is implemented as a piecewise function through a scalar field t i . The flux term is kept unchanged : Numerically and mathematically, there is no difference in flux computation between both methods, thus it is not suprising that transient-state fluxes are correctly evaluated, as seen in Figure 13. The low relative deviation observed in flux computation is due to the boundary issue applied to the normal vector dl i . The back-EMF can be calculated using some numerical derivation method.

ϕ i =

Conclusion & Future Work

This paper presents the first step of an original method which accurately models magnetic devices by avoiding iterative meshing step. Results are very promising as CPU time is reduced in 3D keeping low relative deviation. CPU time saving is performed by avoiding meshing step, and also by a better initialization of iterative resolution. Optimization, movement modeling and transient-state simulation are very efficient and give same results as the usual FE method. Nevertheless, the use of a fixed mesh grid leads to boundaries issues, therefore to some imprecision in very specific cases. Several advices and solutions are given to overcome this problem. The real interest lies in a kind of mixed model using the usual method for static parts and the projection method for moving and deforming parts. In other words, the purpose of future works is to a develop model in which static parts of the system are discretized using the usual method and moving parts are projected on a fine structured mesh, allowing time saving in iterative processes. In conclusion, this paper represents the first step of a promising method which can be used in several domains of finite element method.

Figure 1 :

 1 Figure 1: Principle of the usual method applied on a conventional PMSM

Ω

Figure 3 :Figure 4 :

 34 Figure 3: Nodal linear interpolation of physical properties implying boundary deviation with projection method

Figure 5 :

 5 Figure 5: Nodal linear interpolation of physical properties implying boundary deviation with projection method using non physical line to control mesh element size

Figure 6 :

 6 Figure 6: 2D flux density map for usual method (left) and proposed method (right)

Figure 9 :

 9 Figure 9: Torque comparison for 3D models

Figure 11 :

 11 Figure 11: Robustness of proposed method for 260,000 Dofs

Ωa

  Figure 13: Transient-state magnetic flux computation comparison.

Table 1 :

 1 Time saving

		Usual method Proposed method
	Linear problem
	Meshing time	1	ε * 1/N b position
	Resolution time	1	1
	Non-Linear problem
	Meshing time	1	ε * 1/N b position
	Resolution time	1	1/2

Table 2 :

 2 Geometrical variation

	Variable parameters	Variation range
	Slot opening (%)	[15 ; 85]% of slot pitch
	Magnet opening (%)	[15 ; 85]% of pole pitch
	Armature radius (mm)	[72...108]
	Yoke radius (mm)	[140 ... 172]
	(without magnet)

Table 3 :

 3 Summary of optimization processes

		Problem I	Problem II	Problem III
	Objective	∆ Γ

Table 4

 4

	: Results of optimization processes
			Problem I
		Usual Model Proposed Model
	∆ Γ	1p.u.	1.01p.u.
	CPU Time	52h	38h
	Iterations	115	145
			Problem II
		Usual Model Proposed Model
	∆ Γ	1p.u.	0.99p.u.
	CPU Time	10h	10h
	Iterations	20	40
			Problem III
		Usual Model Proposed Model
	∆ Γ	1p.u.	0.98p.u.
	CPU Time	7.5h	3.8h
	Iterations	15	15


Georg Knebel

Motoi Kimata

Michal Vališka

Fuminori Honda

Dexin Li

Daniel Braithwaite

Gérard Lapertot

William Knafo

Alexandre Pourret

Yoshiki J Sato

Yusei Shimizu

Takumi Kihara

Jean-Pascal Brison

Jacques Flouquet

Dai Aoki

Anisotropy of the Upper Critical Field in the Heavy-Fermion Superconductor UTe 2 under Pressure

We studied the anisotropy of the superconducting upper critical field H c2 in the heavy-fermion superconductor UTe 2 under hydrostatic pressure by magnetoresistivity measurements. In agreement with previous experiments we confirm that superconductivity disappears near a critical pressure p c ≈ 1.5 GPa, and a magnetically ordered state appears. The unusual H c2 (T ) at low temperatures for H a suggests that the multiple superconducting phases which appear under pressure have quite different H c2 . For a field applied along the hard magnetization b axis H c2 (0) is glued to the metamagnetic transition H m which is suppressed near p c . The suppression of H m with pressure follows the decrease of temperature T max χ , at the maximum in the susceptibility along b. The strong reinforcement of H c2 at ambient pressure for H b above 16 T is rapidly suppressed under pressure due to the increase of T sc and the decrease of H m . The change in the hierarchy of the anisotropy of H c2 (0) on approaching p c points out that the c axis becomes the hard magnetization axis.

In many strongly correlated electron systems unconventional superconductivity (SC) appears close to a quantum phase transition where a long range ordered phase is suppressed by tuning a control parameter of the system such as pressure, doping, or charge carrier number. [START_REF] Flouquet | Progress in Low Temperature Physics[END_REF][START_REF] Löhneysen | [END_REF] It is believed that the enhancement of the magnetic and electronic fluctuations are the glue for the superconducting pairing. This has been most impressively shown for ferromagnetic superconductors such as URhGe and UCoGe, where the pairing strength itself can be tuned by the magnetic field. [3][4][5] A magnetic field applied along the easy magnetization axis of these orthorhombic systems lowers the superconducting pairing, while a magnetic field applied along the hard magnetization axis enhances the superconducting pairing strength as the field drives the system to a collapse of the ferromagnetism. This enhancement of the pairing strength for H b results in the reentrance of SC in the field range of 9 T -12 T in URhGe 6) and an enhancement of SC in UCoGe. 7) The strong internal exchange field in these ferromagnetic superconductors and the extremely high ratio of the upper critical field H c2 compared to the superconducting transition temperature T sc imposes a non-unitary spin triplet state with equal spin pairing (ESP) which is a superconducting state very rarely realized in bulk materials.

Recently the superconducting state below T sc = 1.6 K of the paramagnetic heavy fermion compound UTe 2 has also been proposed to be a spin triplet superconductor. 8,9) Evidence for this is obtained from the very large and strongly anisotropic H c2 , which is 6 T for the field applied along the easy magnetization axis a, 11 T for the c-axis but is extremely field-enhanced up to 35 T for field along the hard magnetization b-axis (exceeding by far the Pauli limitation for a singlet superconductor). Above 35 T, where a metamagnetic transition with a huge jump of the magnetization M occurs, 10,11) SC is abruptly suppressed. An additional singularity of UTe 2 is that under pressure multiple superconducting phases occur. 12) T sc is initially linearly suppressed with pressure, but for p > 0.3 GPa two specific heat anomalies are observed and the upper superconducting anomaly increases up to 3 K while the lower continues decreasing in temperature. In that experiment SC is suppressed near 1.7 GPa and a new magnetically ordered phase appears. The increase of T sc by a factor of 2 compared to the the ambient pressure value and the collapse of the superconducting regime near 1.7 GPa has been confirmed by resistivity experiments. [12][13][14] In the present work we concentrate on the anisotropy of H c2 under pressure. First we present the pressure dependence of the susceptibility. We show that (i) H c2 for H a is unusally enhanced at low temperature, (ii) for H b the superconducting phase is suppressed at a metamagnetic transition, which decreases in field by increasing pressure, and (iii) the upper critical field for H c crosses that for H b for p > 1.2 GPa. This could be related to a change in the magneto-cristalline anisotropy and we speculate that the c axis becomes the hardest magnetization axis at high pressure.

Single crystals of UTe 2 have been grown by the chemical vapor transport method as described previously, in Grenoble and in Oarai. 9,15) High pressure susceptibility measurements have been performed in a Quantum design MPMS using a specially designed Cu-Be piston cylinder cell. For the resistivity experiments under pressure in Grenoble, a several millimeter long needle-like single crystal has been cut in three pieces and mounted in a piston cylinder cell. So field could be applied on different parts of the same crystal along the a, b, and c axes. The current was injected along the field direction for H a while for the other directions the current is perpendicular to the field. The measurements have been performed in a Quantum design PPMS (maximal field 9 T) and a dilution refrigerator (H max = 13 T). In parallel, high pressure experiments have been performed in Oarai using an Oxford top-loading dilution refrigerator (H max = 15 T). Results from experiments in Grenoble and Oarai are very similar, except that the samples used in Oarai have a slightly higher T sc under pressure. In addition we performed, at selected pressures, measurements up to 30 T in the high field laboratory in Sendai using 3 He and 4 He cryostats.

In Fig. 1(a) we show the magnetic susceptibility as M/H for a field of 1 T applied along the b axis as a function of temperature for different pressures up to 0.9 GPa. At zero pressure, the susceptibility has a maximum at T max χ ≈ 31.5 K, slightly lower than outside the pressure cell. At zero pressure T max χ is linked to the metamagnetic transition at H m ≈ 35 T. 10,16) It shows that the same energy scale is responsible for the formation of a correlated electronic regime in zero magnetic field and pushes the system under magnetic field to a metamagnetic transition. 17) In UTe 2 a huge jump of the magnetization ∆M = 0.6 µ B is observed at H m at p = 0, while the susceptibility ∂M/∂H is almost unchanged below and above H m . The maximum of the susceptibility T max χ decreases under pressure, and at 0.9 GPa we find T max χ ≈ 20 K. The absolute value of M/H at low temperature is inversely proportional to T max χ and increases under pressure [see Fig. 1(b)]. We have also added the pressure dependence of the metamagnetic transition field H m detected by magnetoresistivity. Importantly it follows the pressure dependence of T max χ . A rough extrapolation yields H m → 0 near 1.5 GPa, i.e. near the pressure where SC is replaced by a magnetically ordered state at p c . [12][13][14] Figure 2 displays H c2 as a function of temperature along the a, b and c axes of UTe 2 . T sc (H) or H c2 (T ) have been determined from temperature sweeps at constant field, or field sweeps at constant temperature using the criterion ρ = 0 (see Supplemental Material for raw data). 18) The inset in Fig. 2(a) shows the pressure-temperature phase diagram. We find good agreement with the previous reports. 12,13) No complete superconducting transition has been found at 1.48 GPa. Thus, compared to Ref.12, the critical pressure p c up to which SC can be observed is slightly lower. We also do not observe the initial decrease of T sc , 12) as our first pressure is already 0.3 GPa. H c2 (T ) for H a is shown in Fig. 2(a). The ambient pressure data are on a different sample. At zero temperature H c2 (0) ≈ 6 T. At the same field a Lifshitz transition of the Fermi surface is observed. 20) Increasing the pressure, T sc increases with a maximum at p ≈ 1 GPa. Most spectacularly we find a strong enhancement of H c2 below 0.6 K at 0.5 GPa and below 0.75 K at 0.86 GPa marked by a kink in H c2 (T ). In some organic superconductors or iron-pnictides a FFLO state is observed at high field. [21][22][23] Here this is certainly not the case, as a FFLO state induces at most a positive curvature of H c2 , not such a kink, leading to a maximum increase of H c2 (0) of around 6% for the pure paramagnetic limit of three dimensional superconductors, not 25% as observed here.It is most likely due to multiple superconducting phases with different order parameters, as in the phase diagrams of UPt 3 , or Th-doped UBe 13 . [24][25][26] Thermodynamic measurements are required to reveal the possible appearance of extra superconducting phases below the superconducting boundary detected here. 27) The strong curvature of H c2 (T ) under pressure at 0.5 GPa and 0.86 GPa points to the presence of a strong Pauli paramagnetic effect. Even at p = 0, H c2 (T ) determined specific heat measurements points to a paramagnetic limitation. 28) At 1.4 GPa H c2 has clearly a reentrant behavior at low field. Increasing the magnetic field, H c2 increases from 1.93 K at H = 0 to 2.1 K at 3 T. It shows that close to the critical pressure, field is enhancing (or restoring) SC in this direction. At lower temperature, we do not see any enhancement of SC at 1.4 GPa in difference to the lower pressures.

In Fig. 2(b) we show H c2 (T ) for H b. The known remarkable feature at p = 0 is the enhancement of SC from 17 T up to 35 T. 19) In the Grenoble experiment, the maximum field has been 13 T, such that we could not follow the field enhancement. However, we clearly see the change of behavior of H c2 (T ) as a function of pressure near the maximum of T sc at p ≈ 1 GPa. While at 0.86GPa and 1.07 GPa H c2 (T ) is almost linear up to the highest field, we observe, for p 1.21 GPa, a marked downward curvature of H c2 near T sc . This strong curvature of H c2 (T ) observed on approaching p c , suggests either a strong Pauli paramagnetic limit, or a strong field dependence of the pairing. As we will see below another underlying effect is the suppression of the metamagnetic transition. Moreover, at 1.32 GPa, we can detect a clear hysteresis of 0.5 T between field sweep up and field sweep down, for temperatures below 1.2 K (roughly 0.5T sc ) [closed and open circles in Fig. 2(b) and Supplemental Material], indicating a first order nature of the superconducting transition.

In Fig. 2(c) we plot H c2 (T ) for H c for different pressures. Strikingly, we see upward curvatures very close to T c for all pressures. It is most pronounced at 0.5 GPa. Such an initial upward curvature may occur in resistivity experiments due to filamentary SC or to multiband effects, with dominant pairing for light carriers. Similarly, in Ref. 12 or in UCoGe 29) such behaviour is reported from resistivity measurements. On the contrary, H c2 determined by ac calorimetry for H c has a linear temperature dependence at T sc . 12) Here, we observe that for fields above 2 T, H c2 (T ) varies linearly with temperature and the slope is strongly enhanced with pressure. In the inset of Fig. 2(c) we compare H c2 for the b axis and c axis at 1.21 GPa. While for the b axis H c2 has the usual downward curvature near T sc (0), it has an upward curvature for the c axis and is linear only above 5 T. Neglecting the initial curvature, the slope for H c is about -90 T/K at 1.21 GPa, comparable to that determined from ac calorimetry close to the critical pressure. 12) This indicates an extreme re-inforcement of SC along the c axis close to p c . 30) Closer insights on the superconducting properties can be obtained by analysing in more details H c2 (T ) close to T sc . In order to extract as precisely as possible the initial slope H ′ c2 = (dH c2 /dT ) T =T sc even, when there is a strong positive curvature, we made a (weak-coupling) fit of the data taking into account both the orbital and a possible Pauli limitation [see Fig. S9 in Supplemental Material]. In Fig. 3(a) we plot this determination of the initial slope H ′ c2 normalised by T sc as a function of pressure for the three crystal directions. In a normal superconductor, i.e. when the pairing strength does not change significantly with the applied field, it is a good measure of the average Fermi velocity v F in directions perpendicular to the applied field, as H ′ c2 ∝ -T sc /v 2 F . So, it is roughly proportional to the square of the corresponding effective mass (m ⋆ ) 2 . In Fig. 3(a) we plot -H ′ c2 /T sc as function of pressure. While up to 1 GPa the H ′ c2 /T sc is almost constant, it increases strongly above 1 GPa. For the b axis we find an increase by a factor 9 between 1 GPa and 1.4 GPa. From Ref. 12 we know that H ′ c2 increases for the c axis until SC is suppressed. Our data here are in good agreement with those determined by ac-calorimetry, and the increase up to p c for H c is by a factor 6. For H a we plot the data only up to 0.86 GPa, but at 1.4 GPa the initial slope has already changed sign. Thus, H ′ c2 is strongly reinforced on approaching p c for all three directions.

Of course, the initial slope is not a direct measure of the effective mass m ⋆ due to the possible field enhancement of the strong-coupling parameter λ. 4,8,19) On the other hand, such a large increase of H ′ c2 for H b certainly gives the right trend for the pressure increase of m ⋆ . The pressure dependence of H ′ c2 compares qualitatively with the pressure dependence of the A coefficient of the T 2 term of the resistivity. At low pressure the resistivity follows well a ρ(T ) = ρ 0 + AT 2 Fermi-liquid temperature dependence, and with increasing pressure the coefficient increases from A = 2.7µΩ cm/K 2 to 17.5 µΩ cm/K 2 at 0.86 GPa by a factor 5 at zero field, see Fig. 3(b). We have also plotted the A coefficient for a field of 13 T which shows that for the b and c axes it is already enhanced, while for the a axis A is almost constant. However, above 1 GPa the resisitivity is almost linear in the normal state above the superconducting transition indicating the closeness to some quantum critical point (see also Supplemental Material).

In the inset of Fig. 3(a) we have plotted the electronic g factor determined from fitting also the curvature of H c2 (T ), which is related to the Pauli limitation 31) Of course, these values are only indicative, because neither strong coupling effects nor a field dependence of the pairing have been taken into account. However, it is clear that for H a and above 1 GPa for H b, the Pauli limitation of H c2 is not negligible. The fact that paramagnetic limitation appears along two perpendicular directions at least puts strong constraints on the possible superconducting order parameters: for example, it is not possible to have a simple real d-vector fixed in a given direction as it would lead to a paramagnetic limitation only along this direction (and ESP states in the perpendicular directions). It could be that we have a complex d-vector (non-unitary state) ?, 28) with no component along the c axis, or an A 1u order parameter (like in the B phase of superfluid 3 He), with a paramagnetic limitation present along the three directions.

H P c2 (0) = √ 2∆ gµ B .
Finally, in Fig. 4(a) we show the phase diagram for H b at 1 GPa. We clearly see that SC is suppressed continuously with field, but is suddenly cut off when H c2 is of the same order as the metamagnetic transition field H m . At 1 GPa we find H m ≈ 20 T, in very good agreement with the maximum in the temperature dependence of the susceptibility (T max χ ≈ 20 K) [see Fig. 1(b)]. The first order nature of the transition at H m is clearly manifested in the jump of the residual resistivity ρ 0 at H m shown in Fig. 4(b) from about 50 µΩcm to 90 µΩcm. The jump ∆ρ 0 ≈ 40µΩcm is three times smaller than at p = 0. However, there are also strong magnetic fluctuations developing at H m connected to a huge increase of the A coefficient. As shown in Fig. 4(c) the relative field dependence of A(H) at 1 GPa scales perfectly with the A(H) at p = 0 [see Fig. 4(d)]. The disappearance of SC right above H m is in contrast to the rather symmetric field dependence of A(H). In Fig. 4(e) we have summarized the H, T phase diagrams for H b for different pressures. We want to stress that the reinforcement of H c2 observed at zero pressure above 16 T occurs near 8 T at p = 0.57 GPa and has completely disappeared at 1 GPa. If it is due to a change in the superconducting order parameter, 32) only the high field phase will survive close to p c . 14) The signature of the metamagnetic transition gradually fades out between 1 GPa and 1.4 GPa, where only a kink in the resistivity is observed, that could be a signature of H m . This explains why no clear cut off of SC at H m appears close to p c .

Metamagnetism in heavy fermion compounds leads generally to a drastic change of the nature of the magnetic correlations with a reconstruction of the Fermi surface. For example in ferromagnetic superconductors this is well established in URhGe for field applied along the b axis. 33) Similar phenomena have been detected on crossing the antiferromagnetic to paramagnetic transition of CeRh 2 Si 2 34) or in UPd 2 Al 3 . 35) A highly studied case is the CeRu 2 Si 2 series (see e.g. [START_REF] Flouquet | Progress in Low Temperature Physics[END_REF]36) ); again, for pure CeRu 2 Si 2 , the sharp crossover at H m from a nearly antiferromagnetic phase to a polarized paramagnetic phase, is associated with a drop of antiferromagnetic correlations above H m and a drastic Fermi surface reconstruction.

A simple "rule of thumb" is that ferromagnetic or antiferromagnetic correlations drop above H m and mainly only local fluctuations survive. Thus in UTe 2 if pairing is indeed due to the FM intersite correlations, it should be strongly suppressed above H m : it is clear for H b that it is strongly suppressed. Furthermore even at H = 0 above p c ≈ 1.5 GPa, it is still not known whether UTe 2 shows ferromagnetic or antiferromagnetic order (there is no clear indication for ferromagnetism in the transport measurements) The anisotropy of H c2 close to p c is reversed compared to that at p = 0, and we find

H c c2 (0) > H a c2 (0) ∼ H b c2 (0)
for the c, a, and b axes, respectively. As the anisotropy of the H c2 (0) in most heavy fermion superconductors (UPt 3 might be the only exception) follows that of the susceptibility, it is likely that the c axis becomes the hard axis near the critical pressure.

Our study of H c2 under pressure for the field along the a, b, and c axes reveals unique features of the superconducting phase of UTe 2 , which clearly demands further thermodynamic experiments. Magnetic field and pressure induced changes of the ferromagnetic interactions as well as possible drastic Fermi surface reconstructions at H m are important ingredients. The duality between the local and itinerant character of the 5 f electrons is important to understand the change in the magnetocrystalline energy leading to a suppression of the metamagnetic field under pressure. SC also evolves strongly in UTe 2 . For H a we have shown that H c2 is unusually enhanced at low temperatures suggesting multiple superconducting phases. 27) The upper critical field for the b axis is determined by the mutual balance of the increase of T sc and the decrease of H m which cuts off SC at the metamagnetic critical field H m . Close to the collapse of SC near 1.5 GPa H c2 (0) is highest along the c axis.

Supplemental Material

In this Supplemental Material we show additional data to those shown in the main article, further data will be made available on special demands.

Resistivity under high pressure

Figure S5 shows the temperature dependence of the resistivity of UT 2 at zero field for different pressures. Up to a pressure of 1.32 GPa zero resisitivty is observed below the superconducting transition temperature T sc . For higher pressures the resistivity at lowest temperature is finite, i.e. no complete superconducting transition is observed. T onset indicates the onset of the superconducting transition. At 1.48 GP the transition is not complete and extremely broad. This indicates that the slope of dT sc /d p is very steep. In addition, already at 1.32 GPa we observe a kink in the resistivity above the onset of the superconducting transition. This is even more pronounced at 1.48 GPa and may be linked to the magnetically ordered state which is observed for higher pressures. This may indicate that the transition from the a superconducting ground state to the magnetically ordered state is first order. We want to emphesize that at 2.13 GPa, the magnetic transition is not at all like expected for a simple ferromagnetic or antiferromagnetic transition and deserves detailed investigatios in future. The shape of the resistivity curves is different than those published in Ref. 13, mainly due to the different direction of the current, which is along the a axis in the present experiment.

H (T) 2.5 K 1.4 K 1.75 K 2.2 K 1.9 K UTe 2 H || a p= 0.86 GPa H c2 1.3 K 1.2 K 1.1 K 0.9 K 0.8 K 0.7 K 0.675 K 0.65 K 0.6 K 0.55 K 0.5 K 0.4 K 0.3 K 0.2 K 0.1 K Fig. S6
. Magnetoresistance for magnetic field applied along the a direction at p = 0.86 GPa. The arrow indicate the criterion ρ = 0 used to determine H c2 .

2. Resistivity at p = 0.86 GPa Figure S6 shows the magnetoresistance ρ(H) as a function of field applied along the a axis of the orthorhombic UTe 2 at p = 0.86 GPa for field along the a axis. The upper critical field H c2 has been defined by the ρ = 0 criterion. In Fig. S7 we show the temperature dependence of the resistivity for the three crystallographic directions for fields between 0 T and 13 T, the current is always injected along the a axis. We checked carefully that there is no dependence on the criterion used to determine H c2 on the overall shape of the upper critical field for H a as shown in Fig. S8. Most remarkably is the strong enhancement of H c2 below 700 mK for H a. These strong increase has also been verified on different samples grown in Oarai and an increase due to some inhomogeneities in the samples can be excluded.

In Fig. S9 we show the analysis of the temperature dependence of the resisitivity for p = 0.86 GPa for fields up to 13 T. We want to stress that the magnetoresistiance in the normal state only for H b axis is positiv, while for the a and c axis it is negativ. The temperature dependence of the resisitivity has been fitted by a power law, ρ(T ) = ρ 0 + AT n in the temperature range between T sc < T < 6 K. The field dependence of the exponent n and the coefficient of the temperature dependent term is shown in Fig. S9. We find that for the resistivity exponent n ≈ 1.7. It is also possible to force the temeprature dependence to a T 2 dependence, however, the fit is significantly less good over this temperature range. In Fig. S9 we show the field dependence of the A n coefficient obtained from the power law. As expected from the magnetoresistance in the normal state, A n is increasing for field along the b direction, while it is decreasing for the a and c direction. We also see, that A n (H) for field along the a axis shows some anomaly around 7 T, which may be connected to the Lifshitz transition observed at zero pressure near 6 T. 20)

Resistivity at p = 1.21 GPa

In Fig. S10 we show the temperature dependence of the resistivity at p = 1.21 GPa for H b (upper panel) and H c (lower panel) for different fields up to 9 T. The magnetoresistance in the normal state is positive for field along the b, and negative for field along the c axis. The temperature dependence is almost linear. Fitting with a power law ρ(T ) = ρ 0 + AT n gives an exponent n = 1.1 for both directions in the normal state. It is obvious that the magnetic field suppresses T sc stronger for field applied along the b axis. The onset of superconductivity is almost field independent in the case of H c, in excellent agreement with the results from ac calorimetry from Ref. 12.

Hysteresis of magnetoresistivity at 1.32 GPa

Near to the critical pressure we observed clear hysteresis in the magnetoresistance as shown in the Fig. S11 for different temperatures. At low temperature the hysteresis is almost 0.5 T. It vanishes near 1.2 K which corresponds to ≈ 0.5T sc .

H -T phase diagram at 1 GPa

Here we show magnetoresistivity data for p = 1 GPa. These data have been used to drav the phase diagram shown in Fig. 4(a) of the main text. The shown data are exemplary for the H -T phase diagrams shown in Fig. 4(e) of the main text.

In Fig. S12(a) we show the magnetoresistivity for different temperatures between 0.5 K and 6 K measured in the high field facility in Sendai. At 0.5 K the resistivity is zero up to H ≈ 20 T. The superconducting transition at 0.5 K is rather broad and only at 25 T the normal state is reached. However, increasing hte temperature the transition get's sharper while the field of zero resistance does not vary up to almost 2 K. For higher temperatures the field of zero resistance decreases, as shown in Fig. S12(b) in the phase diagram. From the magnetoresistivity in the normal state we can identify the metamagnetic trnsition field by the almost steplike increase of the resistivity. As shown for T = 6 K, above the critical end point the step-like increase disappears and a rather broad maximum defined the cross-over temperature which is connected to the temperature of the maximum of the susceptibility, which is here at 1 GPa near 20 K. The important observation is that the metamagnetic transition cuts off the superconductivity defined by the zero resistivity criterium which is nearest to the bulk transition. In Fig.

Analyse of the upper critical field

In Fig. 3 of the main paper we present the pressure dependence of the initial slope H ′ c2 = (dH c2 /dT ) T =T sc normalised by T sc as a function of pressure. In a normal superconductor this allows an estimation of the average Fermi velocity and this of the effective mass of the electrons forming Cooper pairs.In order to get a rather good estimation of the slope, we fitted H c2 near T sc taking into account the Pauli limitation and the orbital limitation. It is the orbital limitation which determines the inital slope while the estimation of the Pauli limitation allows for the description of the curvature of H c2 near T c . Strong coupling effects are neglected. In Fig. S13 we show the upper critical field for various pressures and the fit of the upper critical field from which we determined the slope and also the g factor shown in Fig. 3(a) of the main article. From this fitting it is obvious that for pressures below 1 GPa for field along the b axis, no Pauli limitation is necessary to reproduce the data. However, for higher pressures, a strong curvature occurs. Clearly we see that for H a the Pauli limit can not be neglected. Of course we know that this only a very rough estimation of the behavior and a microscopic model taking into account the correct order parameter, the possible field dependence of the superconducting pairing and the interplay with the underlying metamagnetic transition has to be taken into account.

Fig. 1 .

 1 Fig. 1. (Color online) (a) Temperature dependence of M/H for H = 1 T applied along the b axis of UTe 2 at different pressures. The small upturn below 4 K is due to the background contribution from the pressure cell. (b) Pressure dependence of the temperature of the maximum of the susceptibility T max χ (left red scale, red circles), M/H extrapolated to T = 0 (blue squares, right scale) and H m (left scale, black circles).

Fig. 2 .

 2 Fig. 2. (Color online) Upper critical field H c2 for a magnetic field applied along different directions: (a) H a, (b) H b and (c) H c axis. Open circles are taken from Ref.19. Full circles are from measurements in Grenoble, crosses give H c2 measured at 1.4 GPa in Oarai. The difference in T sc at 1.4 GPa at zero field is probably due to the slight pressure inhomogeneity in the cell. For H a we find a strong enhancement of H c2 at low temperatures indicating the possibility of multi-superconducting phases. The insert in (a) shows T sc as a function of pressure. Dashed lines in (c) are an extrapolation of H c2 from high fields. The insert in (c) shows the crossing of H c2 for H b and c at 1.2 GPa. We used the criterion ρ = 0 to determine H c2 in all cases.

Fig. 3 .

 3 Fig. 3. (Color online) (a) Pressure dependence of the initial slope of H c2 at T sc normalized by T sc for the three crystallographic axes. Inset: Pressure dependence of the g factor indicating the importance of the paramagentic limitation of H c2 . (b) Pressure dependence of the A coefficient of the resistivity at H = 0 (squares). We have also added the variation of the A coefficent at 13 T for H a, b, and c axes.

Fig. 4 .

 4 Fig. 4. (Color online) (a) Magnetic and superconducting phase diagram for p = 1 GPa for H b. SC is abruptly suppressed for H > H m . (b) Field dependence of the residual resistivity ρ 0 (H), which increases suddenly by a factor 2 at H m . (c) Field dependence of the A coefficient, which is strongly enhanced at H m . (d) Comparison of the A coefficient at 1 GPa (blue) to that at p = 0 (red squares). (e) H, T phase diagrams (H b) for different pressures, full symbols correspond to H c2 , open squares to H m or a crossover above the critical end point (CEP), which is indicated by the arrows. Dashed and solid lines are guide to the eye.

Fig. S5 .

 S5 Fig. S5. Resistivity as a function of temperature of UTe 2 at zero magnetic field for various pressures.

  Fig. S7. Temperature dependence of the resistivity for fields between 0 T and 13 T at p = 0.86 GPa for the H a (upper panel), H b (middle panel), H c (lower panel).

Fig. S8 .

 S8 Fig. S8. Upper critical field H c2 at p = 0.86 GPa for the three principal crystallographic directions. Open circles corresponds to the mid-point of the transition for H a.

Fig. S9 .

 S9 Fig.S9. (Upper panel) Field dependence of the resistivity exponent n of the temperature dependence of the resistivity, which has been fitted by a power lar ρ = ρ 0 + AT n . At the pressure p = 0.86 Gpa we find small deviations from the Fermi liquid n = 2. (lower panel) Field dependence of the coefficient A of the temperature dependent term of the resistivity. Remakably, A is decreasing for H a and c, while it is increasing for field along the b axis.

Fig. S10 .

 S10 Fig. S10. (Upper panel) Temperature dependence of the resistivity at 1.21 GPa for H b (upper panel) and H c (lower panel).

  Fig. S11. (Left panel) Magnetoresistance at p = 1.32 GPa for different temperatures. Solid lines are field sweep up, dashed lines corresponds to field down. (Right panel) Zoom on the upper cricital field curve indicating the hysteresis between field sweeps up (solid points) and down (down).

  in Fig.S12(b) we have plotted the superconducting boundary by both, the zero resistivity and the midpoint of the transition. Clearly, bulk superconductivity seems cut of by the metamagnetic transition.

Fig

  Fig. S12. (a) Magnetoresistivity at p = 1 GPa for different temperature. (b) H -T phase diagram for p = 1 GPa. We see that the superconducting transition defined by the ρ = 0 criterium coincides with the metamagnetic transition.

  Fig. S13. (a) Magnetoresistivity at p = 1 GPa for different temperature. (b) H -T phase diagram for p = 1 GPa. We see that the superconducting transition defined by the ρ = 0 criterium coincides with the metamagnetic transition.

Acknowledgment We thank K. Machida, Y. Yanase, V. P. Mineev and Z. Zhitomirsky for fruitful discussions. We acknowledge the financial support of the Cross-Disciplinary Program on Instrumentation and Detection of CEA, the French Alternative Energies and Atomic Energy Commission, KAKENHI (JP15H05882, JP15H05884, JP15K21732,