Péter Forgács

Árpád Lukács

Charge screening in the Abelian Higgs model

In the Abelian Higgs model electric (and magnetic) fields of external charges (and currents) are screened by the scalar field. In this contribution, complementing recent investigations of Ishihara and Ogawa, we present a detailed investigation of charge screening using a perturbative approach with the charge strength as an expansion parameter. It is shown how perfect global and remarkably good local screening can be derived from Gauss' theorem, and the asymptotic form of the fields far from the sources. The perturbative results are shown to compare favourably to the numerical ones.

Historically, the idea for a renormalisable model containing massive gauge fields (of which the Abelian Higgs model is arguably the simplest one) came from the theory of shielding of electromagnetism in a superconductor [START_REF] Anderson | Plasmons, gauge invariance, and mass[END_REF][START_REF] Higss | Broken symmetries and the masses of gauge bosons[END_REF].

The analogy between how the Abelian Higgs model converts a long-range interaction (mediated by a massless gauge field) into a short-range one (mediated by a massive vector field) and how electromagnetism is shielded by mobile charges is well-known. Motivated by the recent investigation of charge screening in the Abelian Higgs model in Ref. [START_REF] Ishihara | Perfect charge screening of extended sources in an Abelian-Higgs model[END_REF], and two subsequent papers, Refs. [START_REF] Ishihara | Homogeneous balls in a spontaneously broken U(1) gauge theory[END_REF][START_REF] Ishihara | Charge screened nontopological solitons in a spontaneously broken U(1) gauge theory[END_REF], where the results are applied, we present a simple, systematic perturbative framework to study this phenomenon analytically, complementing those studies. One is quickly led to believing that all this must have long been done, since, e.g., screening in scalar electrodynamics has been considered in Refs. [START_REF] Mandula | Total charge screening[END_REF][START_REF] Adler | Relaxation methods for gauge field equilibrium equations[END_REF][START_REF] Bawin | Charge screening in classical scalar electrodynamics[END_REF], and in Yang-Mills theories in Refs. [START_REF] Sikivie | Screening solutions to classical Yang-Mills theory[END_REF][START_REF] Campbell | Total source charge and charge screening in Yang-Mills theories[END_REF]. To the best of our knowledge, however, prior to Ref. [START_REF] Ishihara | Perfect charge screening of extended sources in an Abelian-Higgs model[END_REF] there has been no detailed investigation published about this "schoolbook" case.

This note is intended as a (pedagogical) complement of Ref. [START_REF] Ishihara | Perfect charge screening of extended sources in an Abelian-Higgs model[END_REF] where the physics of screening in the Abelian Higgs model has been studied numerically and also analytically by various approximations for both point-like and for extended charge distributions. Our note is tuned to a perturbative study of the physics of charge screening, considering the external charge strength as a "small" parameter and compute corrections to the usual massive Green's function. In this approach it is easy to investigate local charge screening, while global screening is shown to follow from Gauss' theorem [START_REF] Forgács | On non-topological solitons in Abelian gauge theories coupled to U(1)×U(1) symmetric scalar fields[END_REF]. The perturbative results are shown to agree remarkably well with the numerical ones.

The structure of the paper is as follows: in Sec. 1 we introduce the notations used in the model, then in Sec. 2 choose a gauge and a suitable set of variables. Sec. 3 contains the field equations for static solutions, it is exhibited how perfect screening follows from the asymptotic behaviour of the fields and Gauss' theorem and the setup of a series solution. In Sec. [START_REF] Ishihara | Homogeneous balls in a spontaneously broken U(1) gauge theory[END_REF] we consider the field of a point source, using both perturbation theory and numerical methods. In Sec. 5, we calculate the field of a Gaussian charge distribution, in Sec. [START_REF] Mandula | Total charge screening[END_REF], that of a homogeneously charged sphere, and then conclude. Appendix A contains the derivation of the Green function in the static case. Appendix B is devoted to present the Green function of the time-dependent massive Klein-Gordon equation in a closed form in coordinate space representation, and demonstrate its consistency with the Green function of the static equation. In Appendices C resp. D the interaction potential between two point sources in a Klein-Gordon resp. a Proca field are calculated (using the method of Refs. [START_REF] Landau | Course of theoretical physics 2: The classical theory of fields[END_REF][START_REF] Speight | Static intervortex forces[END_REF]).

The model considered

The Abelian Higgs model contains the following fields: a U(1) gauge field A µ , whose mass is to be generated by the Higgs field, a complex scalar field φ, which is subject to self-interaction due to a potential, whose minimum is obtained at a non-zero value of φ. The model is specified by the action integral,

S = d 4 xL , L = - 1 4 F µν F µν + (D µ φ) * D µ φ -V (φ * , φ) -A µ j µ ext , (1)
where there is an implicit summation over repeated Greek indices from 0 to 4, indices are moved up and down with the help of the Minkowski metric, g µν = diag(1, -1, -1, -1) (sign convention of Ref. [START_REF] Landau | Course of theoretical physics 2: The classical theory of fields[END_REF]), Greek indices run from 0 to 3, and Latin (spatial) ones from 1 to 3. In this paper, we shall use units where the unit of length and time agree (velocity of light is unity). The field strength tensor F µν = ∂ µ A ν -∂ ν A µ is antisymmetric, containing the electric and magnetic fields as

F 0i = E i , F ij = ε ijk B k ,
where ǫ ijk is the fully antisymmetric unit tensor in flat 3d space, ε 123 = 1. The gauge covariant derivative of the scalar field is given as D µ φ = (∂ µ -ieA µ )φ, and the potential is

V (φ, φ * ) = λ 4 (φ * φ -η 2 ) 2 ,
where e is the electric charge and η a constant. The vector field j µ ext shall always be a fixed (external) current density, assumed to be conserved, ∂ µ j µ ext = 0. The induced current,

j µ φ = ie [φ(D µ φ) * -φ * D µ φ] , (2)
is also the Noether current corresponding to the symmetry φ → e ieα φ, α a real parameter (U(1) phase symmetry), therefore it is conserved, ∂ µ j µ φ = 0. The Euler-Lagrange equations corresponding to the action (1) are the Maxwell equations for the gauge field A µ , with the source being the sum of the external and the induced currents,

∂ ν F µν = -j µ , j µ = j µ φ + j µ ext , (3)
and the scalar field equation

D µ D µ φ = - ∂V ∂φ * = - λ 2 (φ * φ -η 2 )φ . ( 4
)
2 Gauge choice and identification of degrees of freedom

The action (1) is invariant to transformations of the kind

φ → φ ′ = e ieξ φ , A µ → A ′ µ = A µ + ∂ µ ξ , (5)
where ξ is an arbitrary smooth function of space and time. In the Abelian Higgs model, a customary gauge fixing is the unitary gauge, where the complex scalar field φ becomes real. Parametrising φ as

φ = η + χ √ 2 e iθ , χ = √ 2(|φ| -η) , θ = arg φ , (6)
the unitary gauge corresponds to the choice of a gauge function ξ = -θ/e. Expressing the Lagrangian L with the new variable χ one obtains:

L = - 1 4 F µν F µν + 1 2 m 2 A A µ A µ + 1 2 ∂ µ χ∂ µ χ- λ 16 (χ 2 +2χv) 2 +e 2 vA µ A µ χ+ e 2 2 A µ A µ χ 2 -A µ j µ ext , (7)
where

m 2 A = 2e 2 η 2 , v = √ 2η.
The physical degrees of freedom of the theory are now manifest: A µ has been turned into a massive (Proca) vector field, χ is a massive, self-interacting scalar, and there is a somewhat unusual coupling between the two.

The field equations are

∂ ν F µν -m 2 A A µ = -j ′ µ φ -j µ ext , (8)
(known as the Proca equation), where the induced current is

j ′ µ φ = -2e 2 vχA µ -e 2 χ 2 A µ , (9)
and

∂ µ ∂ µ χ = e 2 A µ A µ (χ + v) - λ 4 (2v 2 + 3vχ + χ 2 )χ . ( 10
)
Taking the divergence of the Proca equation ( 8) and using the antisymmetry of the field strength tensor

∂ µ ∂ ν F µν = 0, one obtains m 2 A ∂ µ A µ = ∂ µ j µ . (11)
Now assuming that the current is conserved, the resulting equation, ∂ µ A µ = 0 agrees formally with the Lorentz gauge condition in electrodynamics. Putting it another way: in the Abelian Higgs model, the unitary gauge in which the scalar field is real, implies also a Lorentz gauge condition on the vector potential. The Proca equation ( 8) is therefore equivalent to

∂ ν ∂ ν A µ + m 2 A A µ = j ′ µ φ + j µ ext , ∂ µ A µ = 0 , (12)
i.e., the massive Klein-Gordon equation with an auxiliary condition (formally, the Lorentz gauge condition). Finally, let us consider the splitting of the equations into three dimensional parts. The Proca electric, E i , resp. magnetic, B i , fields satisfy the Maxwell-type equations

E i = Ȧi -∂ i A 0 , ∂ i E i + m 2 A A 0 = ̺ φ + ̺ ext , (13)
B i = ε ijk ∂ j A k , ε ijk ∂ j B k + m 2 A A i = -Ėi + j ′ φ,i + j ext,i , (14)
where we have introduced the notation ̺ φ = j ′ φ,0 and ̺ ext = j ext,0 . The other two group of Maxwell equations are of the usual form:

ε ijk ∂ j E k = Ḃi , ∂ i B i = 0 , (15)
making it possible to express E i , resp. B i , in terms of the potentials A 0 and A i . The energy density of a field configuration can be written as

E = 1 2 E 2 + B 2 + m 2 A 2 A 2 0 + A 2 + 1 2 χ2 + ∇χ 2 + e 2 2 A 2 0 + A 2 χ 2 + 2vχ + V . ( 16
)
The model may be rescaled by the following change of variables: x µ → x µ /(eη), ρ → ηρ, A µ → ηA µ . The resulting Lagrangian assumes the same form as the original one [Eq. ( 7)], with an overall 1/e 2 factor, and the following replacements: e → e s = 1 in the covariant derivatives, λ → β in the potential, where β = λ/e 2 . In the new units, the vector mass becomes m As = √ 2 and the scalar one m ss = √ β, and the expectation value v s = √ 2. For the integrated quantities, the replacement is Q → (1/e)Q s and E → (η/e)E s , where the rescaled quantities are given by the same formulae as their unscaled counterparts, with the appropriately replaced parameters. In what follows, we shall use these units, and drop the index "s".

Static solutions

Let us now consider time-independent external sources, j µ ext (t, x i ) = j µ (x i ) and seek timeindependent solutions of the Proca equations ( 12) which now become :

-(∇ 2 -m 2 A )A 0 = ̺ φ + ̺ ext , -(∇ 2 -m 2 A )A i = j ′ φ,i + j ext,i , (17)
where ∇ 2 = ∂ i ∂ i . Eq. ( 10) for the scalar field is written as

(∇ 2 -m 2 s )χ = β 4 (3v + χ)χ 2 -e 2 (A 2 0 -A 2 i )(v + χ) , (18)
where m 2 s = βv 2 /2. As the equations for the temporal and the spatial components of Eq. ( 17) agree, it is sufficient to consider one (say, the temporal one) without loss of generality.

For later reference, let us also consider spherically symmetric solutions, where A 0 = A 0 (r) and χ = χ(r). The resulting radial equations are

1 r 2 (r 2 A ′ 0 ) ′ = m 2 A A 0 + 2e 2 A 0 vχ + e 2 A 0 χ 2 -̺ ext , 1 r 2 (r 2 χ ′ ) ′ = m 2 s χ + β 4 (3v + χ)χ 2 -e 2 A 2 0 (v + χ) . (19)
In Eq. ( 19), a prime on radial functions denotes d/dr.

Asymptotic solution

Let us first consider Eq. ( 17) for the case of a localised source, far from the source. Let us introduce spherical coordinates r, ϑ, ϕ, in which the Laplacian ∇ 2 may be written as

∇ 2 = ∂ 2 ∂r 2 + 2 r ∂ ∂r + 1 r 2 ∇ 2 Ω , (20)
where ∇ 2 Ω denotes the angular part of the Laplacian. Considering the equation to O(1/r), the angular part may be neglected. In this order, therefore, the angle dependence of field is arbitrary, and one may seek a solution in the form of an exponential times a series in 1/r. The exponent is determined to be ±m A r, with the negative sign corresponding to a finite energy solution, and the leading power of r is -1,

A 0 = e -m A r r f (ϑ, ϕ) + O(r -2 e -m A r ) . (21)
At this order, the angle dependence is arbitrary; it is, in fact, determined by the source.

Most importantly, the scalar potential A 0 tends to zero exponentially (it is in fact of order O(e -m A r /r)) for r → ∞. Let us now consider Gauss' theorem:

Q = d 3 r(-m 2 A A 0 + ̺ φ + ̺ ext ) = d 3 r∂ i E i = lim r→∞ r 2 π -π dϕ π 0 dϑ sin ϑE r = 0 , (22)
where E r = -∂A 0 /∂r.

The most important consequence of Eq. ( 22) is perfect shielding: the total induced charge (the density of which is -m 2 A A 0 + ̺ φ ) exactly cancels the external charge globally.

Perturbative solution

Let us now consider the case when the sources are "weak",i.e.

ρ ext = ǫ̺ (1) ext , (23)
one can expand the fields in series of ǫ, as

A 0 = ǫA (1) 0 + ǫ 2 A (2) 0 + . . . , χ = ǫ 2 χ (2) + . . . , (24)
[as the external charge is absent from the field equation [START_REF] Campbell | Total source charge and charge screening in Yang-Mills theories[END_REF], there is no linear contribution to the scalar field χ].

In each order, the equations assume the form

(∇ 2 -m 2 s )χ (k) = -ξ (k) , (∇ 2 -m 2 A )A (k) 0 = -σ (k) 0 , (25)
where the source terms ξ (k) , σ (k) 0 are determined by the lower order solutions (up to order (k -1)).

ξ (1) = 0 , σ

0 = ρ (1) (1)
ext , ξ (2) = e 2 vA (1) µ A (1)µ , σ

0 = -2e 2 vχ (1) A (1) 0 , (2)
and the induced charge density, in each order, is given as

-m 2 A A (k) 0 + σ (k) 0 . (27)
It follows from Eqs. ( 19) that χ (1) = 0, and therefore ̺

(1) φ = 0 too. In this order, all of the induced charge comes from the term -m 2 A A

(1) 0 . The solution is obtained in each order with the help of the Yukawa Green function, G i = 1/(4πr) exp(-m i r), i = s, A, in the form

A (k) 0 (x i ) = d 3 x ′ G A (x i -x ′ i )σ (k) 0 (x ′ i ) , G A (x) = 1 4π|x| exp(-m A |x|) , χ (k) (x i ) = d 3 x ′ G s (x i -x ′ i )ξ (k) (x ′ i ) , G s (x) = 1 4π|x| exp(-m s |x|) . (28)
The charge in subsequent orders is calculated as

Q (1) = Q (1) ext + Q (1) ind , Q (k) = Q (k) A + Q (k) φ , (29)
where

Q (1) ext = d 3 xρ (1) ext , Q (k) φ = d 3 xσ (k) 0 , (30)
and

Q (k) A = -d 3 xm 2 A A (k) = -m 2 A d 3 xd 3 x ′ G A (x i -x ′ i )σ (k) 0 (x ′ i ) = -Q (k) φ , (31)
and similarly Q

(1)

A = -Q (1)
ext , where the last equality in Eq. ( 31) is obtained by performing the integration over x first.

The procedure for obtaining the full solution in the form of the series in Eq. ( 24) order by order, each order having a source from the lower ones is referred to as a dressing procedure. Eq. (31) shows, that there is perfect global charge screening, Q

(k) φ + Q (k) A = 0 in each order.
The Yukawa Green function used in Eq. ( 28) also tells us something about local charge cancellation. The contribution of a point charge, if linear approximation suffices, is screened within a sphere of radius O(1/m A ). If the external sources do not change significantly on this scale, the external charge density is cancelled by the induced charge to a very good accuracy.

Another length scale stems from the scalar, 1/m s . In the terminology of superconductivity, the length scale 1/m A is termed penetration depth and 1/m s is the correlation length. This determines the nature of the interaction between point particles at large distance: if the penetration length is larger, the vector interaction dominates, and the interaction is repulsive, and, on the contrary, if the correlation length is larger, scalar interaction dominates, and the interaction is attractive. The same holds for the interaction between flux tubes in superconductors, and, therefore, this determines the magnetic properties of superconductors; see, e.g., Ref. [START_REF] Tinkham | Introduction to superconductivity[END_REF].

Spherical symmetry

Let us also consider spherical symmetry for the series solution, i.e., when ρ exr = ρ ext (r). In this case, the A 0 and χ are also spherically symmetric. Introducing the shorthand notation y 25) can be compactly written as

(k) i i = A, s for A (k) 0 resp. χ (k) Eqs. (
1 r 2 r 2 y (k) i ′ ′ -m 2 i y (k) i = -h (k) i , (32)
where h

(k) s = ξ (k) and h (k) A = σ (k) 0 .
Eqs. (32) are second order inhomogeneous linear differential equations, which can be solved using the two linearly independent solutions, y i± = e ±m i r /r, of the respective homogeneous equations, and their Wronskian

W i = y i+ y ′ i--y ′ i+ y i-= -2m i /r 2 , (i = A, s) as y (k) i (r) = e -m i r r r r 0 dx e m i x h (k) i (x) W i (x) - e m i r r r ∞ dx e -m i x h (k) i (x) W i (x) , (33)
where the integration constants (limits) are chosen to ensure the boundary conditions. The condition that both A 0 and χ tend to 0 for r → ∞ is clearly implemented in Eqs. (33).

The boundary conditions at r = 0 come from regularity in the sense that terms ∝ 1/r at r → 0 be absent (otherwise they would yield unwanted Dirac-delta sources in ∇ 2 A 0 or ∇ 2 χ), this fixes the value of the integration constant (r 0 ).

Point source

A point source is defined by a Dirac delta as the external charge density,

̺ ext (x) = qδ (3) (x) . (34)
In this case, we shall use q as the expansion parameter, ǫ = q.

Assuming spherical symmetry the leading term is

A (1) 0 (r) = G A (r) = e -m A r 4πr . ( 35
)
The energy of the point particle solutions obviously diverges since for r → 0, G ∼ 1/(4πr), the electromagnetic field contribution to the energy density, E 2 /2 ∼ 1/(32π 2 r 4 ), thus E = 4π Er 2 dr, which is divergent.

The second order contribution to the scalar field, χ (2) (r) is obtained from Eq. (33) yielding

χ (2) (r) = - e 2 v 2(4π) 2 m s r e -msr Ei[(m s -2m A )r] -log |m s -2m A | m s + 2m A -e msr Ei[-(m s + 2m A )r] , Ei(x) = - x -∞ dt e t t . (36)
wheredenotes the principal value integral, see Refs. [START_REF] Abramowitz | Handbook of mathematical functions[END_REF][START_REF]The NIST handbook of mathematical functions[END_REF]. In deriving χ (2) (r) we have also used the expansion Ei(x) ∼ -γ + ln(|x|) + x for x → 0. From (36) it follows that for r → 0 χ (2) (r) ∝ log(r).

In the special case m s = 2m A , the solution in (36) is replaced by

χ (2) (r) = - e 2 v 2(4π) 2 m s r
e -msr ln(r/r 0 ) -e msr Ei[-2m s r] , r 0 = e γ 2m s .

(37)

Interaction energy

The interaction energy between two point charges q 1 , q 2 placed at r 1 resp. r 2 at a separation r = r 1 -r 2 in both (massive) electrodynamics and in Klein-Gordon theory is calculated as ±1 times the product of the field of one charge at the position of the second one multiplied by the value second charge (see Refs. [START_REF] Landau | Course of theoretical physics 2: The classical theory of fields[END_REF][START_REF] Iwanenko | Klassische Feldtheorie[END_REF]and Appendices C and D).

In the case m s > m A , it is the field due to A

(1) 0 which dominates for large r, yielding

V II (r) = q 1 q 2 4πr e -m A r , (38)
whereas for m A > m s , it is the scalar field χ (2) which dominates. In this case, for r → ∞, the leading contribution is given by

χ (2) (r) ∼ e 2 v 2(4π) 2 m s r e -msr log |m s -2m A | m s + 2m A . ( 39
)
The field χ (2) in Eq. ( 39) is the field of a point source in a Klein-Gordon field with strength e 2 v/2/(4π)/m s log |m s -2m A |/(m s + 2m A ).

To obtain the interaction potential between q 1 and q 2 , the field of charge q 2 is approximated as emanating from a point source determined by the near field limit (this approximation has been used by Ref. [START_REF] Speight | Static intervortex forces[END_REF] for the interaction energy of vortices in superconductors). That is ξ (2) is replaced by q (2) s δ 3 ( r 2 ), where q (2) s = -d 3 xξ (2) = -e 2 v/(8πm A ). q (2) s multiplies the field of the charge q 1 , yielding

V I (r) = e 4 v 2 q 2 1 q 2 2 4(4π) 3 m s m A log 2m A -m s 2m A + m s e -msr r . (40)
The indices "I" and "II" to distinguish m A > m s and m s > m A (β < or > √ 2) were choosen in accord with the terminology of superconductivity [START_REF] Tinkham | Introduction to superconductivity[END_REF]. In both cases, the force is exponentially decreasing with the separation; however, in the case of a type I setting, like charges attract (the logarithm is negative). Formulae (38) and (40) are valid for large separations r.

Numerical calculations

The Dirac delta source does not appear in the radial equations, which are only defined for r > 0. Instead, it manifests itself in the boundary conditions [START_REF] Adler | Theory of static quark forces[END_REF][START_REF] Adler | Relaxation methods for gauge field equilibrium equations[END_REF][START_REF] Ishihara | Perfect charge screening of extended sources in an Abelian-Higgs model[END_REF], as r → 0, A 0 ∼ q/(4πr), i.e., very close to the source, it is unshielded. For a series solution, the leading power for a scalar field is determined from considering the coefficient of the lowest power, yielding

r → 0 : A 0 ∼ q 4πr , χ ∼ χ 0 r γ , (41)
with γ = -1/2 + 1/4 -κ 2 (for κ < 1/2), where κ = eq/(4π) and

r → 0 : A 0 ∼ q 4πr , χ ∼ χ 0 1 √ r cos(γ 1 log r + δ) , (42)
where γ 1 = κ 2 -1/4 for κ > 1. In eqs. ( 41) and (42), χ 0 is a constant determined from the numerical solution of the radial equations. For r → 0, all fields approach their vacuum value,

r → ∞ : A 0 → 0 , χ → 0 . (43)
We have calculated numerical solutions using the Colnew package [START_REF] Ascher | A new basis implementation for a mixed order boundary value ODE solver[END_REF][START_REF] Ascher | Numerical solution of boundary value problems for ordinary differential equations[END_REF]. An example numerical solution is diplayed in Fig. 1a, and the corresponding charge in Fig. 1b. For the given value of the charge q = 0.4 and self-interaction β = 2.0, the leading order series agrees extremely well with the exact result (within line width of Fig. 1a).

In Fig. 2 we have plotted the radius R c , defined as

Q φ (R c ) = -0.9q , (44)
as a function of q and β. We have found, that R c is a decreasing function of q and an increasing one of β.

Gaussian source

Another interesting example is a Gaussian source, with source radius R s , Figure 2: Point source: the dependence of the radius of screening with a precision of 0.9 on (a) the charge q and (b) the strength of self-interaction β (for q = 0.4).

ρ ext = q (2π) 3/2 R 3 s e -r 2
Again, we shall use ǫ = q as the expansion parameter. The integration constant in the solution given in Eq. ( 33) is chosen to ensure regularity for r → 0 (i.e. to cancel the 1/r terms), yielding

A (1) 0 = 1 4πr e -m A r e m 2 A R 2 s 2 2 Erfc m A R 2 s -r √ 2R s -e 2m A r Erfc m A R 2 s + r √ 2R s , ( 46
)
where Erfc(x) = 1 -

x 0 dx ′ exp(-x 2 ) = 1 -Erf(x)
is the complementary error function (see Ref. [START_REF] Gradshteyn | Table of integrals, series[END_REF]).

The expression for the potential A

(1) 0

of the Gaussian source in Eq. ( 46) starts with the potential G of the point source, see Eq. (57). It is instructive to introduce an effective charge as

Q eff = q lim r→∞ A (1) 0 G = qe m 2 A R 2 s 2 ≈ q 1 + m 2 A R 2 s 2 + . . . , (47)
where we have used the fact, that Erf(x) ∼ 1 + exp(-x 2 )/( √ πx) for x → ∞ (7.1.23 and 7.12.1 in Refs. [START_REF] Abramowitz | Handbook of mathematical functions[END_REF][START_REF]The NIST handbook of mathematical functions[END_REF], resp.), i.e., the leading contribution comes from the constant 2 term in the asymptotics of the first Erfc function in Eq. ( 46). According to Eq. ( 47), at a large distance, the potential due to an extended source is stronger than that of a point charge of the same magnitude. This is physically plausible, as for massless electrodynamics, they agree, and here, some part of the source is less strongly screened. We have also calculated numerical solutions for the extended source distribution given by Eq. (45), using the Colnew package [START_REF] Ascher | A new basis implementation for a mixed order boundary value ODE solver[END_REF][START_REF] Ascher | Numerical solution of boundary value problems for ordinary differential equations[END_REF]. A typical solution is shown in Fig. 3. On the figures, the half-width of the source R s and the radius R c where the screening reaches an accuracy of 0.9 are also shown. Again an excellent agreement is found between the first order analytical solution (46) and the numerical one depicted in Fig. 3a, the difference between the two curves being smaller than the line width. We have also used the numerical solutions to obtain R c as a function of q, and have found that it is monotonically decreasing (Fig. 4a) and of β (Fig. 4b), of which it is weakly increasing. Our results are in perfect agreement with the numerical results of Ref. [START_REF] Ishihara | Perfect charge screening of extended sources in an Abelian-Higgs model[END_REF]. In the case of closed form formulae, our analyses complement each other. Ref. [START_REF] Ishihara | Perfect charge screening of extended sources in an Abelian-Higgs model[END_REF] consider the case of a small source, R s ≪ 1/m A , and in this case present an expansion of A 0 in r, and a large source, R s ≫ 1/m A an in this case show that charge cancellation is local, and

A 0 ≈ ρ ext /m 2
A in the leading order. These approximations are found to be of good agreement with numerical results. Our analysis, on the other hand, provides the expansion of the solution in the external charge q, and is also found to be in excellent agreement with numerical results.

Homogeneous sphere

Let us also consider a homogeneous sphere as the source, ρ ext = ρ 0 Θ(R s -r), where Θ(x) is the Heaviside function, Θ(x < 0) = 0, Θ(x > 0) = 1. In this case, obtaining the scalar potential Figure 4: Gaussian source: the dependence of the radius of screening with a precision of 0.9 on (a) the charge q and (b) the strength of self-interaction β (for q = 0.4). In both cases, R s = 1/10/ √ 2.

using formula (33) yields A

= ρ 0 /q 2m 3 A r e -m A r -e m A r< (1 -m A r < ) + e -m Rs (1 + m A R s ) + e m A r Θ(R s -r) e -m A r (1 + m A r) -e -m A Rs (1 + m A R s ) , (1) 0 (r)
where r < = min(r, R s ). Again, we have used the total charge, q = 4πR 3 s ρ 0 /3 as the expansion parameter. The agreement with the numerical result is again near perfect, within the line width of Fig. 5a. The perturbative result, Eq. (48) allows one to introduce an effective charge again, by comparing the asymptotic form of the solution to Q eff /(4πr)e -m A r , which yields

Q h,eff = 2πρ 0 m 3 A e -m A Rs (1 + m A R s ) -e m A Rs (1 -m A Rs) ≈ q 1 + m 2 A R 2 s 10 + . . . , (49)
again, the effective charge is slightly larger than the total external charge. The behavior of the screening radius R c on the parameters is similar as before: it is weakly decreasing with q and increasing with β.

In Ref. [START_REF] Ishihara | Perfect charge screening of extended sources in an Abelian-Higgs model[END_REF], Eq. ( 48) has been derived with a matching procedure. In addition, they consider the case of a lage sphere, in which case the fields inside the sphere are constant, and can be found by neglecting derivatives in the radial equations [START_REF] Adler | Theory of static quark forces[END_REF]. The latter approximation is also used in Ref. [START_REF] Ishihara | Perfect charge screening of extended sources in an Abelian-Higgs model[END_REF] for a thick walled spherical source. Ref. [START_REF] Ishihara | Perfect charge screening of extended sources in an Abelian-Higgs model[END_REF] finds good agreement between the approximate formulae and numerical results. Figure 6: Homogeneous sphere source: the dependence of the radius of screening with a precision of 0.9 on (a) the charge q and (b) the strength of self-interaction β (for q = 0.1). In both cases, R s = 1/ √ 2.

Conclusions

The screening of time-independent external charges in the Abelian Higgs model has been considered, from multiple aspects. On one hand, using the asymptotic form of the gauge field far from the charges, it has been shown, from Gauss' law, that the charge screening is exact, the charges induced in the fields of the Abelian Higgs model globally cancel the external charge exactly. The same result has been demonstrated using perturbation theory considering the charge as the perturbation parameter. Local cancellation occurs if the external charge distribution does not change significantly on the scale of 1/m A , where m A is the gauge boson mass. We have obtained simple explicit formulae for the fields in the first nontrivial order in perturbation theory and compared them to numerical results, and found very good agreement. We note that our results are in perfect agreement and complementing those of Ref. [START_REF] Ishihara | Perfect charge screening of extended sources in an Abelian-Higgs model[END_REF].

is calculated, again, with the help of a Fourier transformations, based on the discussion in Ref. [START_REF] Iwanenko | Klassische Feldtheorie[END_REF],

G 4 (x, t) = dk 0 d 3 k (2π) 4 e -i(k 0 t-kx) G4 (k 0 , k) . (59)
Using the Fourier representation of the Dirac delta, G4 = -1/(k 2 0k 2 -m 2 ) is obtained, and, using spherical coordinates for k, with the 3rd axis aligned along r, the form

G 4 = dk 0 dkk 2 d(cos ϑ)dϕ (2π) 4 -1 k 2 0 -k 2 -m 2 e -ik 0 t e ikr cos ϑ (60)
is obtained. The integral over k 0 is performed using the theorem of residues. For t > 0, the contour may be closed with a large semicircle on the lower, whereas for t < 0, in the upper halfplane. Causality demands, that the poles at k 0 = ± √ k 2 + m 2 be shifted slightly downwards, so that G(x, t) = 0 for t < 0 results. For t > 0, the integration along the contour is clockwise, pole contibutions are multiplied by 2π/i;

G 4 = Θ(t) dkk 2 d(cos ϑ)dϕ (2π) 3 e ikr cos ϑ sin( √ k 2 + m 2 t) √ k 2 + m 2 , (61)
where Θ(t) is the Heaviside theta function, Θ(t > 0) = 1 and Θ(t < 0) = 0. The integrals over ϕ yields a factor of 2π and the one over d cos ϑ a sine function, after which we transform the product of two sines into two cosines, replace k → -k in the second one, and obtain

G 4 = Θ(t) 4π 2 r ∞ -∞ dkk 2 √ k 2 + m 2 cos √ k 2 + m 2 t -kr = - Θ(t) 4π 2 r ∂ ∂r ∞ -∞ dk √ k 2 + m 2 sin √ k 2 + m 2 t -kr . (62)
The last integral in Eq. ( 62) can be evaluated with a change of variable for t > r and another one for r > t. In both cases, we change the variable k as k = m sinh κ, so that √ k 2 + m 2 = m cosh κ and dk = m cosh κdκ.

In the case t > r, we also transform as t = u cosh τ , r = u sinh τ , u = √ t 2 -r 2 (invariant interval). This way,

I = ∞ -∞ dk √ k 2 + m 2 sin √ k 2 + m 2 t -kr = ∞ -∞ dκ sin[mu cosh(κ -τ )] = πJ 0 (mu) , (63)
where J 0 denotes the regular Bessel function of order 0 (see Refs. [START_REF] Gradshteyn | Table of integrals, series[END_REF][START_REF] Prudnikov | Integrals and series, vol. 1: Elementary functions[END_REF]).

In the case r > t, we transform as t = u sinh τ and r = u cosh τ , obtaining

I = dk √ k 2 + m 2 sin √ k 2 + m 2 t -kr = ∞ -∞ dκ sin[mu sinh(κ -τ )] = 0 , (64)
as the integrand is odd.

The results of the integration can be summarised as

I = Θ(t -s)πJ 0 (mu) . (65)
Upon derivation w.r.t. r, we obtain

G(x, t) = δ(t -r) 4πr - mΘ(t -r) 4π √ t 2 -r 2 J 1 m √ t 2 -r 2 . (66)
In Ref. [START_REF] Iwanenko | Klassische Feldtheorie[END_REF], the calculation is performed for arbitrary (integer) dimension of spacetime.

A good consistency check of our calculation is now comparing the Green function of the time dependent Klein-Gordon equation, and that of ∇ 2 -m 2 obtained in Sec. 3.2. What one needs to show is that

dtG 4 (x, t) = G(x) . (67)
The Dirac delta term yields 1/(4πr). In the second term, due to the Θ function, we need to integrate over t from r to infinity. By a change of variable r = rt ′ , we are led to the integral (Refs. [START_REF] Gradshteyn | Table of integrals, series[END_REF][START_REF] Prudnikov | Integrals and series[END_REF])

∞ 1 J 1 (b √ x 2 -1) √ x 2 -1 dx = 1 b (1 -e -b
) .

(68)

C Interaction potential in the Klein-Gordon model

Let us consider interaction energy in the massive Klein-Gordon model, i.e., a massive real scalar field, coupled to external Dirac-delta sources. The field satisfies the equation (for the static case)

(∇ 2 -m 2 s )φ = σ = s 1 δ (3) (x) + s 2 δ (3) (x -x 0 ) . ( 69
)
Using the linearity of the equation, we split the field as φ = φ 1 + φ 2 , both satisfying

(∇ 2 -m 2 s )φ 1 = σ 1 = s 1 δ (3) (x) , (∇ 2 -m 2 s )φ 2 = σ 2 = s 2 δ (3) (x -x 0 ) . (70)
The energy density of a static Klein-Gordon field is

E KG = 1 2 |∇φ| 2 + 1 2 m 2 s φ 2 + σφ . (71)
The interaction energy is therefore the cross term in the energy of φ 1 + φ 2

E KG,int = ∇φ 1 ∇φ 2 + m 2 s φ 1 φ 2 + σ 1 φ 2 + σ 2 φ 2 , (72)
and the interaction potential between the two particles corresponding to the Dirac delta sources is

V KG = d 3 xE KG,int , (73)
where, in the first term, we apply the identity ∇φ 1 ∇φ 2 = ∇(φ 1 ∇φ 2 ) -φ 1 ∇ 2 φ 2 . Here, upon integration, the first term gives a vanishing contribution,

d 3 x∇(φ 1 ∇φ 2 ) = φ 1 ∇φ 2 d 2 s = 0 , (74)
where the second interation is taken over a large sphere of radius R, in the limit R → ∞, and d 2 s is the surface element vector; the fields vanish exponentially, therefore, the integral vanishes. What remains, cancels φ 1 σ 2 , yielding

V KG = d 3 xφ 2 σ 1 = s 1 φ 2 (0) , (75)
Note, that had although we have used the field equation for φ 2 , the result is symmetric. Had we used the vector identity as ∇φ 1 ∇φ 2 = ∇(∇φ 1 φ 2 ) -∇ 2 φ 1 φ 2 , and the field equation for φ 1 from Eq. (70), we would have arrived at the same result.

D Interaction potential in a Proca field

A static, purely electric Proca field with two point sources, one at zero, satisfy the field equation

(∇ 2 -m 2 A )A 0 = -ρ = -q 1 δ (3) (x) -q 2 δ (3) (x -x 0 ) , (77)
and the energy of the static, purely electric Proca field is

E P = - 1 2 (∇A 0 ) 2 - 1 2 m 2 (A 0 ) 2 + ρA 0 . (78)
We shall apply a similar splitting procedure as in the case of the Klein-Gordon field, A 0 = A 1,0 + A 2,0 ,

(∇ 2 -m 2 A )A 1,0 = -ρ 1 = q 1 δ (3) (x) , (∇ 2 -m 2 A )A 2,0 = -ρ 2 = q 2 δ (3) (x -x 0 ) . ( 79
)
The interaction energy, as before, is

E P,int = -∇A 1,0 ∇A 2,0 -m 2 A 1,0 A 2,0 + ρ 1 A 2,0 + ρ 2 A 1,0 , (80)
and with the same partial integration technique, using now the Green function of -∇ 2 + m 2 , G = 1/(4πr)e -m A r , V P = q 1 q 2 4πr e -m A r (81) is obtained.

Figure 1 :

 1 Figure 1: (a) The profile functions A 0 and χ of a solution for point charge; β = 2.0, q = 0.4. (b) The charge distributions of the same solution. The dashed vertical line shows r = R c .

Figure 3 :

 3 Figure 3: (a) The profile functions A 0 and χ of a solution for the Gaussian charge distribution; β = 2.0, q = 0.4, R s = 0.1/ √ 2. (b) The charge distributions of the same solution. The dotted vertical line shows r = R s and the dashed one r = R c .

Figure 5 :

 5 Figure 5: (a) The profile functions A 0 and χ of a solution for the homogeneous charge distribution; β = 2.0, q = 0.4, R s = 1/ √ 2. (b) The charge distributions of the same solution. The dotted red line shows r = R s and the dashed blue one r = R c .

using

  Eq. (70) the properties of the Dirac delta. Now, using the Green function of the operator ∇ 2 -m 2 2 , G = -1/(4πr) exp(-m s r), we obtainV KG = -s 1 s 2 4πr e -msr , r = |x 0 | .

Acknowledgements Á.L. acknowledges the support of the Spanish Ministerio de Ciencia, Innovación y Universidades (Grant No. PCI2018-092896) and the EU (QuantERA CEBBEC).

A The Yukawa Green function

The linearised equation is

with ψ being A (n) 0

or χ (n) and σ its respective source, and m = m A or m s . The solution of Eq. ( 50) may be constructed from the solution of the equation

in the form

The function G is the Green function of the operator ∇ 2 -m 2 . It may be computed by seeking it in the form of a Fourier integral,

Using the fact, that δ(x) = dk/(2π) exp(ikx), we learn that

We may choose coordinates in k-space in such a way, that its k z axis is along the direction of the vector x, and use spherical coordinates, resulting in

where r = |x|. The integral over ϕ merely cancels one of the 2π factors. The integral over ϑ can be performed by noting that d cos ϑ = -sin ϑdϑ, and flipping the limits, yielding

Considering now complex values of k, and closing the integration contour with a large semicircle on the upper half-plane, its contribution vanishes exponentially, because r > 0. The direction is positive (counter-clockwise), therefore one needs to add up the residues of poles on the upper half-plane, with a prefactor of 2πi. There is one such pole, at k = im, yielding

which is nothing else than the Yukawa potential.

B The retarded Green function of the massive Klein-Gordon equation

In this appendix, the Green function of the massive Klein-Gordon equation, satisfying


Michal Hegedüs

Marcela Achimovičová

Hongjue Hui

Gabin Guélou

Pierric Lemoine

Ismail Fourati

Jean Juraszek

B Malaman

Peter Baláž

Emmanuel Guilmeau
email: emmanuel.guilmeau@ensicaen.fr

Promoted crystallisation and cationic ordering in thermoelectric Cu26V2Sn6S32 colusite by eccentric vibratory ball milling

A pristine colusite Cu26V2Sn6S32 was successfully synthesised on a 100 g scale via mechanochemical reaction in an industrial eccentric vibratory ball mill followed by spark plasma sintering (SPS) at 873 K. The milling of elemental precursors from 1 up to 12 hours was performed and the prepared samples were investigated in details by X-ray powder diffraction, Mössbauer spectroscopy, scanning electron microscopy, and thermoelectric property measurements. The results point to the formation of a high purity and high crystallinity non-exsoluted colusite phase after SPS process (𝑃4 # 3𝑛, a = 10.7614(1) Å) in the case of a 12 h A c c e p t e d M a n u s c r i p t milled sample. In comparison, samples milled for 1-6 h displayed small quantities of binary Cu-S phases and vanadium core-shell inclusions, leading to a V-poor/Sn-rich colusite with higher degree of structural disorder. These samples exhibit lower electrical conductivity and Seebeck coefficient while an increase in the total thermal conductivity is observed. This phenomenon is explained by a higher reactivity and grain size reduction upon prolonged milling and by a weak evolution of the chemical composition from a partly disordered V-poor/Sn-rich colusite phase to a well-ordered stoichiometric Cu26V2Sn6S32 colusite, which leads to a decrease in carrier concentration. For all samples, the calculated PF values, around 0.7-0.8 mW m -1 K -2 at 700 K, are comparable to the values previously achieved for mechanochemically synthesised Cu26V2Sn6S32 in laboratory mills. This approach thus serves as an example of scaling-up possibility for sulphur-based TE materials and supports their future large-scale deployment.

Introduction

Waste heat harvesting through conversion into electricity by thermoelectric (TE) materials is essential for future energy management. [START_REF] Rowe | Thermoelectrics Handbook: Macro to Nano[END_REF][START_REF] Shi | [END_REF][3] The advent of new low-toxicity and cost-effective inorganic TE materials offers opportunities to recover waste heat more efficiently and economically using highly reliable systems that neither possess any moving parts nor produce any noise or vibration. Generally, best thermoelectric performances are obtained in metal tellurides, which are expensive and toxic. [START_REF] Rowe | Thermoelectrics Handbook: Macro to Nano[END_REF]4 To conciliate efficiency with environmental and cost constraints, it is required to replace tellurium (Te) with a less-toxic and more earthabundant homologous element. Compared with Te, sulphur (S) has a lighter mass, smaller atomic radius, and larger electronegativity. The lighter mass should lead to higher phonon frequencies/Debye temperature and hence higher lattice thermal conductivity (κL). Thus, sulphides are naively thought to have a slightly higher thermal conductivity, and a lower power factor, 𝑃𝐹 = 𝑆 ! 𝜌 "# , and hence a lower figure of merit, 𝑍𝑇 = 𝑆 ! 𝑇𝜌 "# 𝜅 "# , where S is the A c c e p t e d M a n u s c r i p t Seebeck coefficient, r the electrical conductivity, κ the thermal conductivity, and T is the absolute temperature. As a consequence, TE sulphides have been somewhat overlooked with respect to tellurides and selenides despite their more significant diversity of crystal structures and compositions. Thus, research focusing on TE sulphides has recently regained momentum. 5 Cu-based ternary and quaternary sulphides with complex structures were especially identified as potential TE materials for mid-temperature range applications. 6 They can exhibit very low thermal conductivity combined with a semiconducting or degenerate semiconducting (metallic-like) behaviour, allowing to reach relatively high ZT values. A non-exhaustive list includes: tetrahedrites Cu12Sb4S13, [7][8][9][10][11][12][13] , colusites Cu26T2M6S32 (T = V, Nb, Ta, Cr, Mo, W; M = Sn, Ge), 14- 20 germanite derivatives Cu22Fe8Ge4S32, 21 stannoidite Cu8Fe3Sn2S12, 22 bornite Cu5FeS4, [23][24][25][26] chalcopyrite CuFeS2, 27 isocubanite CuFe2S3, 28 mohite Cu2SnS3, 29 kesterite Cu2ZnSnS4, 30,31 and Cu4Sn7S16. 32 Among those materials, colusites with the general formula Cu26T2M6S32 (T = V, Nb, Ta, Cr, Mo, W; M = Sn, Ge) can exhibit high PF values close to 2 mW m -1 K -2 , 20, 33 and ZT values close to 1 at 700 K 19,20,34 (See references given in Table 1). Colusites have a complex cubic structure (space group 𝑃4 # 3𝑛, a ≈ 10.6-10.8 Å) with 66 atoms per unit cell, which can be described as a sphalerite-type network of corner-sharing CuS4 and MS4 tetrahedra (Figure 1). These tetrahedra are ordered in such way that the "sphalerite" unit cell is doubled in all 3 directions and the T cation occupies the 2a interstitial position at the corners and centre of this 2aSph×2aSph×2aSph supercell.

The thermoelectric performance can be controlled to a large extent thanks to cationic substitution on different crystallographic sites or by altering the stoichiometry, as shown in Table 1. It was especially shown that the reduced lattice thermal conductivity (kL) in colusites was attributed to several kinds of phonon scatterers, i.e. Cu/Sn antisite defects, interstitial A c c e p t e d M a n u s c r i p t defects, antisite defects, and split Cu site, [34][35][36][37] mainly induced by the processing conditions used. Based on these observations, it is essential to point out that the processing step plays a key role in determining the final thermoelectric properties, 34,36 especially in the case of structurally complex sulphides. Reproducibility and scalability of their synthesis represent another challenge that is yet to be addressed by current research in this field. The mechanochemistry seems to be the most effective, low-cost, and eco-friendly approach towards scalable synthesis of Cu-based sulphides. 38,39 This has been recently demonstrated on ternary and quaternary copper sulphides, such as tetrahedrite Cu12Sb4S13, [40][41][42][START_REF] Weller | [END_REF][44] Cu4Sn7S16, 32 bornite Cu5FeS4, 24,25 colusite Cu26V2Sn6S32, 16,20,34,36,45,46 germanite derivative Cu22Fe8Ge4S32, 21,47,48 stannoidite Cu8Fe3Sn2S12, 22 isocubanite CuFe2S3, 28 mohite Cu2SnS3, 49,50 stannite Cu2FeSnS4, 51 kesterite Cu2ZnSnS4, [52][53][54][55][56] and mawsonite Cu6Fe2SnS8, 57 that can be readily prepared by ball milling in a planetary mill by direct reaction of elements or binary sulphides. The possibility to prepare pristine colusite by mechanochemical synthesis was first reported by Bourgès et al. 16 , but only for small batches of 5g. With the aim to synthesise large quantity of colusite in one step and with high reproducibility, we investigated the scalability of pristine colusite Cu26V2Sn6S32 synthesis using an industrial eccentric vibratory ball mill. The present study reports on the effect of milling time on the sample purity, crystallinity and thermoelectric performances of pristine colusite Cu26V2Sn6S32. Nihon Seiko, Japan), and sulphur (99% CG-Chemikalien, Germany) were used as starting reactants. The amounts of each element were calculated based on the Cu26V2Sn6S32 stoichiometry.

After completion of the milling programs, the resulting pulverised powder was shaped and densified using SPS (FCT HPD 25) at 873 K for 30 minutes (heating and cooling rate of 100 K/min) under a pressure of 64 MPa using graphite dies of 10 mm diameter and a slight overpressure of 30 hPa (Ar), in order to prevent both oxidation and sulphur volatilisation. The final thickness of the pellets was around 8 mm with a geometrical density in excess of 90 % of the crystallographic value.

X-ray powder diffraction (XRPD) patterns of the milled samples were acquired at room temperature using a Bruker D8 Advance Vario 1 two-circle diffractometer in q-2q Bragg-Brentano geometry equipped with a Cu X-ray tube, a Ge (111) monochromator (i.e. l = 1.54056 Å) and a Lynx Eye detector. The XRPD data were analysed by Rietveld refinement using the FullProf and WinPlotr software packages. 59,60 119 Sn Mössbauer spectrometry was performed on polycrystalline powders at room temperature in transmission geometry and constant-acceleration mode. The source was 119m Sn source in CaSnO3 matrix with an activity of ~130 MBq. For the velocity calibration, a 57 Fe Mössbauer spectrum of a a-Fe foil at RT was collected using a 57 Co source mounted at the other side of the spectrometer transducer. Isomer shifts (d) are given with respect to SnO2 at room temperature. Least-squares fits of the experimental data were performed using the histogram method assuming Lorentzian lines.

Observations of fractured cross-sections of densified pellets by scanning electron microscopy (SEM) were performed using a ZEISS Supra 55 field emission scanning electron microscope.

A c c e p t e d M a n u s c r i p t

The electrical resistivity (ρ) and Seebeck coefficient (S) were measured simultaneously on ingots, from 300 K up to 700 K using an ULVAC-ZEM3 device under partial helium pressure. A NETZSCH LFA-457 apparatus was used for measuring the thermal diffusivity under argon flow. The thermal conductivity (κ) was determined as the product of the geometrical density, the thermal diffusivity, and the theoretical heat capacity calculated by the Dulong-Petit approximation. The Wiedemann-Franz law, using a Lorenz number estimated from the relationship L = 1.5 + exp(-|S|/116), 61 was used to calculate the lattice thermal conductivity by subtracting the electronic contribution to the total thermal conductivity (κL = κ -κe). The estimated measurement uncertainties are 6% for the Seebeck coefficient, 8% for the electrical resistivity, 11% for the thermal conductivity, and 16% for the final figure of merit,

ZT. 62

Results and discussion

The XRPD patterns of the elemental powders milled for a period of time ranging from 1 h to 12 h are shown in Figure 2. First, it should be noted that, after only one hour of milling, no diffraction peaks arising from elemental precursors (Cu, V, Sn, S) are observed. This is closely connected to highly negative Gibbs formation energies of binary and ternary sulphides, especially in Cu-S based systems. 63 Rapid formation of Cu2-xS and Cu2SnS3 has been recently reported. 49,50 These nanocrystalline sulphides are formed within the very first minutes of milling in an explosive way. Secondly, the relatively low crystallinity of the powders leads to high full width at half maximum (FWHM) values of diffraction peaks. Consequently, it is not possible to observe the low angle diffraction peaks characteristic of the colusite structure and therefore to confirm formation of the colusite phase during the ball-milling process only. Indeed, the main diffraction peaks observed on these XRPD data can be attributed to a (Cu,V,Sn)S sphalerite-derivative phase such as poorly crystalline colusite Cu26V2Sn6S32, or mohite Cu2SnS3, phases. Finally, diffraction peaks corresponding to a covellite CuS phase are observed on the XRPD patterns, especially for short milling times. In previous studies, 49,50 it was shown that, during ball-mill processes, binary and ternary phases are formed first, followed by more complex sulphides in later stages of milling, annealing or sintering. Hence, the most likely admixture in samples after ball-milling is covellite CuS with either mohite Cu2SnS3 or a cubic (Cu,V,Sn)S sphalerite-derivative. However, due to the significant overlap of Bragg reflections of these phases with colusite, the formation of the latter cannot be excluded, especially at longer milling times. Importantly, traces of vanadium metal and Cu3VS4 sulvanite are detected in the powder samples (arrows) with shorter milling times (1h, 3h, 6h), but are absent in the sample milled during 12h. As discussed below in SEM section, this has a huge impact over the purity of the sintered products, as well as chemical composition and final crystallinity of colusite phase.

A c c e p t e d M a n u s c r i p t sulphide even after 1 hour of milling followed by SPS treatment (Figure 3a). When the mixture is subjected to longer milling periods, the shape of the diffraction peaks related to colusite structure evolves. Indeed, the integral intensity of the superstructure 211, 321, and 332 reflections evolves from 1 to 6 hours of milling and seems stable for higher milling times. This can be explained by a slight evolution of the chemical composition of the colusite phase coupled with an increased (micro)structural ordering with increasing milling time (see SEM section).

This is supported by a decrease in the unit cell parameter of the colusite phase with increasing milling duration, [33][34][35]37 (≈ 5 wt.% and 1 wt.%, respectively, Figure 3a) and 3h+SPS samples (≈ 2 wt.% and 3 wt.%, respectively, Figure S1a), and covellite in 6h+SPS sample (≈ 1 wt.%, Figure S1b). It should be mentioned that XRPD data analysis indicates the contamination of the 6h+SPS sample by SiO2 quartz of aragonite structure (≈ 1 wt.%) typically used for cleaning of the milling equipment, and of the 12h+SPS sample by traces of WC phase (< 1 wt.%) arising from abrasion of the milling medium (Figure 3b). Nevertheless, these results indicate that 12 hours of milling allow to produce large quantity (100 grams) of powder, which, after spark plasma sintering at 873 K, leads to a high purity pristine colusite material (Figure 3b).

A c c e p t e d

M a n u s c r i p t

A c c e p t e d M a n u s c r i p t

In order to better understand the influence of milling times on the structural disorder/crystallinity of colusite, 119 Sn Mössbauer spectroscopy measurements were performed at room temperature. This technique allows to obtain important information on the local environment of tin atoms in solids (electronic structure and chemical bonding), 64 and consequently, indicates indirectly the degree of structural order-disorder in Sn-colusite. 35 119 Sn Mössbauer spectra of selected sintered samples are shown in Figure 4. The spectrum of the 1h+SPS sample (Figure 4a) exhibits a single peak centred near 1.5 mm s -1 , with a shape distinct from a pure Lorentzian line. For the best fit of the experimental data, we use a combination of a singlet and a quadrupolar doublet (green and red curves in Figure 4a, respectively), as for the fitting procedure described by Candolfi et al. 35 Both components have the same refined isomer shift value (d = 1.45 mm s -1 ) which is characteristic of +IV oxidation state of Sn atoms, and almost the same relative spectral area (49% and 51%, respectively). These two contributions indicate the presence of Sn atoms on more than one crystallographic site in the crystal structure, one of them being slightly distorted, as evidenced by the refined value of the quadrupole splitting of the doublet (ΔEQ = 1.02 ± 0.07 mm s -1 ). The same fitting procedure has been applied

for the Mössbauer spectrum of the 6h+SPS sample (Figure 4b), exhibiting a line with narrower shape. The refined value of the quadrupole splitting of the doublet is reduced (ΔEQ = 0.83 ± 0.18 mm s -1 ) and the relative spectral area of the single line is increased (66 %), in agreement with a progressive ordering of the colusite phase with milling time as evidenced by XRPD analyses. This is further corroborated by the results obtained for the 12h+SPS sample for which the Mössbauer spectrum (Figure 4c) can be fitted with only a singlet with refined isomer shift and linewidth values of d = 1.48 ± 0.01 mm s -1 and Γ = 1.10 ± 0.07 mm s -1 , respectively. This indicates that almost all the Sn atoms are located on a high symmetry site such as the 6c crystallographic site, 35 and confirms the high quality of the produced colusite material at the highest milling time. Note that, in the case of 1h+SPS and 12h+SPS samples, a doublet (~2%

A c c e p t e d M a n u s c r i p t of total spectral area, blue curve in Figures 4a and4c) has been added to reproduce the slight shoulder of the peak near the zero-velocity position. This doublet is attributed to traces of SnO2 impurity with a volume fraction too small to be detectable by XRPD. with the XRPD data and the rapid formation of nanocrystalline sulphides in less than 1 h of ball-milling. After SPS, the grain morphology appears isotropic, as expected from the cubic crystal structure of colusite, and the grain size ranges between 1 and 5 µm (Figure 5). It can be observed that the grain size tends to decrease with increasing milling time, which is self-evident in Figure 5d for 12h+SPS sample. Such feature is consistent with the more complete crystallisation in 12h+SPS sample due to the higher reactivity of the powders, as discussed in the structural analysis section above. Then, the formation of larger grains is favoured by reactive sintering in powders subjected to a ball-milling period shorter than 6h.

Back-scattered electron microscopy (BSE) coupled with EDS analysis confirms the complete crystallisation and homogeneity of the colusite phase for the 12h+SPS sample (Figure 6, right).

The BSE-EDS analysis clearly shows the presence of vanadium-rich core-shell structures in 1h+SPS (Figure S3, left), 3h+SPS (Figure S3, right) and 6h+SPS (Figure 6, left) samples, consistent with our previous work, 36 and confirming the established formation mechanism. The latter has been demonstrated to be a rapid reaction between the Cu, Sn and S, forming a vanadium-poor colusite matrix that progressively accommodates the remaining vanadium metal. A core-shell structure of vanadium/sulvanite/V-rich/Sn-poor colusite is formed, surrounded by a V-poor/Sn-rich colusite matrix (Figure S4). It must be noted that the size of the vanadium particles tends to decrease from 1h+SPS to 3h-6h+SPS samples, and disappear in the 12h+SPS sample, in agreement with the XRD patterns presented in colusite matrix. This is also consistent with the cationic elemental analysis from EDS (Table S5), in which the vanadium content within the matrix progressively increases from a V-poor/Snrich colusite to the expected Cu:V:Sn ratio of 26:2:6. Note that the presence of vanadium metal and surrounding layers of Cu3VS4 sulvanite in SPS samples was missed by the XRD analysis because of the scarcity of vanadium and the close structural relationship between colusite and sulvanite, the main reflections of the latter overlapping those of the former. Colusite is also known to stabilise despite significant deviations from the ideal stoichiometry 26:2:6:32. 65 The crystallisation of Sn-rich colusites for short milling times (£ 6h) is also in good agreement with the Mössbauer data (Figure 4), which indicate the presence of Sn atoms on more than one crystallographic site in the crystal structure and the evolution of the unit cell parameter of colusite determined by XRPD refinements (Table 2). This highlights the importance of milling time investigations for the fabrication of chalcogenides to favour complete crystallisation during the subsequent synthesis/sintering step. The temperature dependences of the electrical resistivity and Seebeck coefficient are displayed in Figure 7. All exhibit a p-type degenerate-semiconducting behaviour in agreement with firstprinciple electronic structure calculations and experimental data. 14,34 The magnitudes of both the electrical resistivity and Seebeck coefficient are similar for the 1h+SPS, 3h+SPS and MPa). 16,34 These electrical properties confirm the relatively good purity and crystallinity of the samples even for a ball-milling time as short as 1h, in agreement with the XRPD and SEM analyses. In the 12h+SPS sample, the electrical resistivity and Seebeck coefficient slightly increase over the whole temperature range, consistent with a slight decrease in the carrier concentration. This trend is probably linked to a higher purity and reactivity of the 12h ballmilled powder leading to a complete crystallisation of a more stoichiometric colusite phase during SPS and the absence of the vanadium core-shell structures observed by SEM for short milling times (£ 6h). This hypothesis is strongly supported by the smaller grain size in this sample (after SPS), as shown in Figure 5d, and the BSE-EDS analyses presented in Figure S3-6 and Table S5. As expected, the heightened reactivity of the 12h ball-milled powder further leads to a better crystallinity, as confirmed by Rietveld refinements of XRPD data and 119 Sn Mössbauer spectroscopy in the previous sections. Note that the measured charge carrier concentration, considering the experimental error, does not vary significantly between the samples, remaining between 8.5×10 21 and 9×10 21 cm -3 , due to the slight variation of Seebeck coefficient between 1h, 3h, 6h, and 12h+SPS samples. Overall, the changes in Seebeck coefficient and in the electrical resistivity with milling time compensate each other and PF values, around 0.7-0.8 mWm -1 K -2 at 700 K (Fig. 7e), are retained for all samples. These values are comparable to those achieved in small batches of Cu26V2Sn6S32 colusites. 16,34 Heating and cooling measurements (displayed in Figure S5) of the electrical resistivity and Seebeck coefficient in ZEM3 revealed that the measurements in 1h and 3h+SPS samples are not reversible, due to the presence of a large quantity of secondary phases and "less reacted"

A c c e p t e d M a n u s c r i p t colusite powders. In contrary, the reversibility is enhanced in 6h and 12h+SPS samples. Further in-depth investigations are needed to address the long-term chemical stability of the as-prepared materials using industrial eccentric vibratory ball-milling. Especially, the chemical sensitivity of colusite against impurities from grinding media (WC, silica) must be studied in details.

The temperature dependence of the thermal conductivity for all samples is presented in Fig. 7c.

All evidence the same temperature dependence characterised by a decrease in the thermal conductivity from 300 K to 700 K. For 1h+SPS, 3h+SPS and 6h+SPS samples, the thermal conductivities are comparable over the whole temperature range, in agreement with the electrical transport properties. The thermal conductivity drops down significantly in 12h+SPS sample as a result of a significantly reduced electronic contribution. Indeed, the lattice thermal conductivities (Fig. 7d) are found rather comparable regardless of the milling time. This suggests that the decrease in grain size, down to 1 µm in 12h+SPS sample, as well as the vanadium core-shell inclusions, are not sufficient to enhance phonon scattering through an increased density of grain boundaries. The values of the total and lattice thermal conductivities are also in good agreement with the results obtained in small batches of Cu26V2Sn6S32 colusites. 16,34 Finally, the figure-of-merit ZT ranges between 0.3 and 0.35 at 700 K among the investigated samples, which is comparable with values reported for undoped colusites synthesised by laboratory mechanochemical synthesis and densified at relatively low temperature (873 K). 16,34,36 The sintering of such colusite powders at high temperature is now required to promote structural disorder [34][35][36] and to achieve lower lattice thermal conductivity and higher ZT values in large ingots.

Conclusions

The pristine colusite Cu26V2Sn6S32 was successfully synthesised by mechanical alloying of elemental precursors in an eccentric industrial vibratory ball mill under argon atmosphere at 100 grams scale. Ball-milled powders showed quintessential polydisperse character with mostly submicronic particles. The need for prolonged milling before SPS treatment, in order to achieve higher crystallinity and purity materials, is emphasised. This is closely related to the enhanced

Figure 1 .

 1 Figure 1. Representation of the crystal structure of pristine colusite Cu26V2Sn6S32.

A c c e p t e d M a n u s c r i p t

A c c e

  p t e d M a n u s c r i p t

Figure 2 .

 2 Figure 2. XRPD patterns of ball-milled elemental mixtures at milling times ranging from 1 h

Figure 3 .

 3 Figure 3. Rietveld refinements of XRPD patterns of (a) 1h+SPS and (b) 12h+SPS samples.

Figure 4 .

 4 Figure 4. 119 Sn Mössbauer spectra at 295 K of (a) 1h+SPS, (b) 6h+SPS and (c) 12h+SPS

Figure 2 .

 2 For short milling times (£ 6h), the powders are insufficiently pre-reacted during ball-milling and the SPS step does not complete the crystallisation of the colusite phase, leaving artefacts of the formation mechanism in the form of vanadium-rich core-shell structures and off-stoichiometric

A c c e p t e d M a n u s c r i p t

A c c e p t e d M a n u s c r i p tFigure 5 .

 5 Figure 5. Micrographs of the fractured cross-sections of the sintered pellets for different ball-

A c c e p t e d M a n u s c r i p tFigure 6 .

 6 Figure 6. BSE micrographs and corresponding EDS elemental analysis maps for the samples

A c c e

  p t e d M a n u s c r i p t 6h+SPS samples, and range around 30 µV K -1 and 0.4 mW cm at RT, respectively. Those values are in good agreement with the results obtained for small batches (5g) of colusite Cu26V2Sn6S32 mechanochemically synthesised in a laboratory mill (Fritsch Pulverisette 7 Premium Line) of elemental precursors, and densified by SPS using similar conditions (873 K/45 min/64

A c c e p t e d M a n u s c r i p tFigure 7 .

 7 Figure 7. Temperature dependence of the (a) electrical resistivity (r), (b) Seebeck coefficient

A c c e

  p t e d M a n u s c r i p t reactivity of milled powders during SPS treatment. The obtained ZT value of bulk colusites using powders synthesised in an industrial mill is comparable to the values obtained in laboratory mills. Overall, the applied approach can serve as a scalable route for the mass production of ternary and quaternary sulphides for a wide range of applications.

A c c e

  p t e d M a n u s c r i p t Y. Kikuchi, Y. Bouyrie, M. Ohta, K. Suekuni, M. Aihara and T. Takabatake, J. Mater. Chem. A, 2016, 4, 15207-15214. 18 Y. Bouyrie, V. Ohorodniichuk, S. Sassi, P. Masschelein, A. Dauscher, C. Candolfi and B. Lenoir, J. Electron. Mater., 2017, 46, 2684-2690.

A c c e p t e d M a n u s c r i p t

Table 1 .

 1 Review on the influence of the preparation method and chemical composition on the thermoelectric performances (represented by ZT values) of colusites. HP: Hot-Pressing; MA: Mechanical Alloying; SPS: Spark Plasma Sintering.

	Preparation method	Nominal composition	ZT	Reference
	melt + HP	Cu26V2M6S32 (M = Ge, Sn)	0.55/0.73 @663 K	Suekuni et al. 14
	melt + HP	Cu26-xZnxV2Sn6S32 (x=0-3)	0.20-0.55 @660 K	Kim et al. 15
	melt + HP	Cu26-yV2Sn6-zS32 (y=0-2; z =0-1)	0.30-0.62 @660 K	Kim et al. 15
	MA + SPS	Cu26-xZnxV2Sn6S32 (x = 0-2)	0.3-0.4 @700 K	Bourgès et al. 16
	melt + HP	Cu26T2Sn6S32 (T = Nb, Ta)	0.6-0.8 @670 K	Kikuchi et al. 17
	melt + SPS	Cu24M'2V2Ge6S32 (M' = Ni, Co)	0.42/0.5 @690 K	Bouyrie et al. 18
	melt + HP	Cu26T2M6-xS32 (T = Nb, Ta; M = Sn, Ge, x=0, 0.5)	0.6-1.0 @670 K	Bouyrie et al. 19
	melt + HP	Cu24Nb2V2Ge6M'xS32 (M' = Co, Ni, Fe; x < 1)	0.55-0.7 @665 K	Bouyrie et al. 58
	MA + SPS	Cu26V2Sn6S32	0.93 @675 K	Bourgès et al. 34
	MA + SPS	Cu26T2Ge6S32 (T = Cr, Mo, W)	0.65-0.86 @700 K	Pavan Kumar et al. 20
	MA + SPS	Cu26V2Sn6S32	0.70 @673 K	Guélou et al. 36
	Experimental Section
	Mechanochemical solid state syntheses were performed in an industrial eccentric
	vibratory ball mill ESM 656-0.5 ks (Siebtechnik, Germany) working under the following
	conditions: 5 L steel satellite milling chamber attached to the main corpus of the mill, tungsten
	carbide balls with a diameter of 35 mm and total mass of 25 kg, 80% ball filling, amplitude of
	the mill 20 mm, rotational speed of the eccenter 960 min -1 , argon atmosphere, total feed of
	reaction precursors 100 grams per batch. The milling was performed from 1h to12 hours.
	Elemental powders of copper (99% Merck, Germany), vanadium (Alfa Aesar 99.5%), tin (99%

Table S1 -

 S1 respectively from 10.7819(1) Å to 10.7614(1) Å for 1h+SPS and 12h+SPS samples. All structural data and reliability factors are provided in Table2and S4. These results confirm the important role played by both milling times and SPS process on phase crystallisation. Finally, these patterns also evidence extra diffraction peaks of weak

intensity indicating the presence of secondary phases, corresponding to digenite Cu2-xS (𝐹𝑚3 # 𝑚, a ≈ 5.56 Å) and covellite CuS (𝑃6 $ /𝑚𝑚𝑐, a ≈ 3.79 Å, c ≈ 16.36 Å) phases in 1h+SPS

Table 2 .

 2 Structural data and reliability factors from Rietveld refinements of the ball-milled colusite Cu26V2Sn6S32 samples after SPS process.

	Sample	1h + SPS	3h + SPS	6h + SPS	12h + SPS
	a (Å)	10.7819(1) 10.7634(1) 10.7634(1) 10.7614(1)
	V (Å 3 )	1253.39(2) 1252.44(2) 1246.95(1) 1246.25(2)
	RBragg	5.51	6.57	3.85	2.59
	RF	5.14	8.03	5.35	4.32
	Rwp	11.30	15.60	10.30	9.47
	Rexp	8.03	10.57	8.04	7.32
	c 2	1.97	2.18	1.65	1.67

Acknowledgments

The authors (M.A., P.B.) acknowledge support of the Slovak Research and Development Agency (project APVV-0357-18). E.G. and G.G. would like to thank Christelle Bilot and Jérôme Lecourt for technical support and acknowledge the financial support of Institut Carnot ESP (project MASTER), FEDER and Normandy Region.