Lower (negative) bounds on the static electric susceptibility of non-equilibrium cubic crystals

We use a classical, microscopic model of point-like dipolarizable entities (a model that is standard in the case of positive polarizability) and investigate its behavior for simple cubic (sc), body-centered cubic (bcc), and face-centered cubic (fcc) crystals with one entity per primitive cell when the static polarizability of the entities is negative and the mutual interaction between the entities is taken into account. We find that the static electric susceptibility is bounded below due to an instability towards self-polarization but the lower permissible bound is negative definite in each case, i.e., the concept of negative static electric susceptibility remains robust, according to the model, when mutual interactions are taken into account. The usual Clausius-Mossotti relation between the static polarizability and the static electric susceptibility remains valid in the case of negative parameters, but only down to the lower permissible bound; the value of the bound depends on the crystal structure and is always unrelated to the asymptote of the Clausius-Mossotti curve. The lower permissible bounds of the static electric susceptibility are found to be -0.906 for sc, and -1.00 for bcc and fcc. These results confirm that, although the magnitude of the static electric susceptibility does not diverge in the negative case (as it can in the positive case), the magnitudes attainable in the negative case for condensed media may, nevertheless, be many orders of magnitude greater than those predicted previously for inverted vapors and gases. This is a promising result in relation to the development of potential new technologies that exploit the phenomenon.


I. INTRODUCTION
The discovery and development of materials with new or enhanced properties is an important driver of technological advance.In the area of electrostatics and low-frequency electromagnetics, a key material property is static electric susceptibility χ (0) (or, equivalently, static relative permittivity ε (0) , related via ε (0) = χ (0) + 1).The value of χ (0) for a material can be crucial to the utility of the material in certain device applications; for example, its value is of crucial importance in capacitor dielectrics, where a large (positive) value is often desired to facilitate device miniaturization [1,2].Since the first measurements of Cavendish [3] and Faraday [4], it was found that the presence of a material between two conductors always increased the static mutual capacitance of the conductors above the capacitance observed when no material was present, i.e., it became established empirically that ε (0) > 1 and, correspondingly, χ (0) > 0 for all materials [5].Although early researchers kept an open mind regarding whether, as-yet untested, materials may nevertheless be found subsequently to exhibit ε (0) values less than unity (see, e.g., Ref. [6]), a theoretical argument was developed later that appeared to rule out this possibility; in §14 of Ref. [7], Landau et al. conclude that "the permittivity of all bodies exceeds unity, and the dielectric susceptibility ... is therefore positive".(Note that Landau et al. are referring here to the static relative permittivity and electric susceptibility).It has since been noted, however, that Landau et al.'s argument assumes that the bodies in question are in thermodynamic equilibrium, and, therefore, the conclusion does not necessarily hold for media that are not in thermodynamic equilibrium [8][9][10][11][12][13].Sanders discussed tentatively the possibility of a χ (0) < 0 state in media with inverted populations of energy levels produced by means similar to those used in maser and laser applications [9], and Chiao et al. predicted unequivocally a χ (0) < 0 state in such systems [10][11][12][13][14]. Apart from the unusual property of reducing, rather than increasing, the static mutual capacitance of two conductors if such a material were placed between them, the possibility of negative χ (0) -a qualitatively new parameter range for this important material property-opens up the possibility of new physical effects and technological capabilities.
For example, Sanders [9] and Chiao et al. [10][11][12][13] discussed theoretically how negative static electric susceptibility opens up the possibility of stable electrostatic levitation, which would be analogous in many respects to the magnetostatic levitation that is seen using diamagnetic or superconducting materials [15][16][17][18].In particular, Chiao et al. proposed ex-plicitly that it should be possible to construct purely electrostatic charged particle traps, which would be very different in principle from Paul and Penning traps [10][11][12].However, Chiao et al. also went on to predict theoretically that for a specific, typical case of the types of systems they considered (ammonia gas at a temperature of 180 K and a pressure of 1 Torr with population inversion maintained by a carbon dioxide pump laser), the value of χ (0) is expected to be ≈ −3 × 10 −4 [11,12]; based on this relatively-small magnitude they concluded that condensed media with much larger magnitudes of the negative static electric susceptibility would need to be developed before such levitation effects become readily observable [11,12].This is perhaps one reason why no attempts to observe experimentally the theoretically-predicted χ (0) < 0 state, or the associated levitation phenomenon, have been reported using the systems considered by Sanders and Chiao et al.
More recently, it was proposed that negative static electric susceptibility may be achieved also in completely different systems: active metamaterials [19].Active metamaterials [20][21][22] utilize an internal source of power and, like the systems considered by Sanders and Chiao et al., are not subject to the equilibrium thermodynamical argument of Landau et al. concerning the restriction on the sign of χ (0) .In Ref. [19], a design concept for active metamaterials with negative static electric susceptibility was proposed and preliminary experimental evidence in support of the general concept was reported.Unlike the systems considered by Sanders and Chiao et al., such metamaterials are readily realized at room temperature and pressure.Further, they constitute a form of condensed matter for which it appears very likely that negative static electric susceptibility values with magnitudes much greater than the systems of Sanders and of Chiao et al. are possible [19].
Given that: (1) the likely existence of negative static electric susceptibility in nonequilibrium materials raises the theoretical possibility of novel technological capabilities [8][9][10][11][12][13], (2) the magnitude of the negative static electric susceptibility that is achievable in non-equilibrium materials is expected to be crucial to the practical realization of such technologies [11,12], and (3) the practical development of condensed (meta)materials with negative static electric susceptibilities of relative large magnitudes is now well under way [19], it appears timely to seek to determine rigorously just how negative can one expect the negative static electric susceptibility to be made in non-equilibrium condensed materials.If Ref. [19] it was claimed-without full theoretical or experimental justification-that negative static electric susceptibility is possible throughout the range −1 < χ (0) < 0; herein we present the theoretical basis on which this claim was made.

A. Scope
We consider the basic, conventional interpretations of χ (0) and ε (0) ; that is, we consider the linear static electric susceptibility (polarization proportional to internal electric field) and linear static relative permittivity (electric displacement proportional to internal electric field) as they pertain to a nonrelativistic, macroscopic, and homogeneous sample of material under the action of a static electric field created by external test electrodes.This interpretation corresponds, for example, to the original experiments of Cavendish [3] and Faraday [4], to the meaning ascribed to the static electric susceptibility and permittivity in standard textbook accounts of the electrostatics of dielectrics (e.g., Chap.II of Ref. [7]), and to the definition of the permittivity according to recent ASTM standards [23].
It should be noted that there are a number of scenarios where negative static electric susceptibility has been discussed or implied in the literature in relation to phenomena that do not correspond to this conventional interpretation.For example, Kirzhnits et al. have shown that static permittivity may be negative in the sense that, if spatial dispersion is taken into account, the longitudinal permittivity at zero frequency but nonzero wave vector may exhibit negative values [24][25][26][27].However, this concerns the scenario where external sources of charge are located within the material itself, and, for the case of static fields and external test electrodes, i.e., the conventional case which corresponds to zero frequency and zero wave vector, Kirzhnits et al. reaffirm the conclusion of Landau et al. that the static permittivity of a material must exceed unity and hence the static electric susceptibility must be positive (if the body is in thermodynamic equilibrium).
We emphasize that we are interested in the static electric susceptibility and permittivity; it is well known that the static magnetic susceptibility may be either positive or negative, regardless of whether the body is in thermodynamic equilibrium or not (see, e.g., Ref. [7], p. 106).We also emphasize that we are interested in the static electric susceptibility and permittivity; it is well known that the sign of the real part of the complex permittivity is subject to no theoretical restriction for nonzero frequencies, again, regardless of whether the body is in thermodynamic equilibrium or not (see, e.g., Ref. [7], p. 274).
Finally, we note that we consider herein only isotropic media such that the static elec-tric susceptibility-which is, in general, a real symmetric second-rank tensor-reduces to a real scalar, which we denote as χ (0) .Thus, 'negative static electric susceptibility' means, straightforwardly, that χ (0) < 0 [28].We use superscript '(0)' to denote explicitly a static quantity.
To study the static electric susceptibility, as interpreted conventionally in this way, we use a simple model of insulators that is standard in the normal case of χ (0) > 0. To analyze the model we also use methods employed previously to a substantial extent in the case χ (0) > 0.
The main novelty of our work is that new results are obtained by applying a generalized version of the method to the unusual parameter range χ (0) < 0.
B. A known lower bound on the value of χ (0) for non-equilibrium materials As noted above, the theoretical argument of Landau et al. ( §14 of Ref. [7]) puts a lower bound of zero on the static electric susceptibility of materials that are in thermodynamic equilibrium.A simple theoretical argument that puts a less stringent lower bound of −1 on the static electric susceptibility-but which applies manifestly to all materials, whether equilibrium or non-equilibrium-is provided by circuit theory, as follows [29].
Consider an empty parallel plate capacitor of capacitance C 0 (> 0) for which the plates are initially disconnected and a charge of ±q 0 is deposited on the plates.Now fill the volume between the capacitor plates with a homogeneous and isotropic insulator with static relative permittivity ε (0) .Suppose that the geometry of the capacitor is such that fringing fields are sufficiently small that they may be reasonably ignored.By definition of ε (0) [23], the capacitance of the filled capacitor is C = ε (0) C 0 .Now let the capacitor plates be connected at time t = 0 across a resistor of resistance R (> 0).Suppose that the value of R is such that the time constant t c = |RC| is large with respect to the time taken for polarization to be established in the material, and the resulting 'quasi-static' behavior of the capacitor is dictated by ε (0) .Elementary circuit theory tells us that the charge on the capacitor as a function of time is q(t) = q 0 exp [−t/ (RC)].For the normal case, ε (0) > 1 and therefore C > 0, hence q(t) decays exponentially, i.e., the capacitor discharges through the resistor.However, for the hypothetical case ε (0) < 0, C would be negative and an exponential, unbounded increase of q(t) would be predicted.Since this indicates a system that is unstable, values ε (0) < 0, i.e., values χ (0) < −1, must be unphysical.We note that the argument does not rely, at any point, on the material being in thermodynamic equilibrium, so it should apply to all materials whether or not they are in thermodynamic equilibrium [30].(viz., all values < −1) in non-equilibrium materials, and it already tells us that one cannot hope to obtain negative values with the large magnitudes readily seen in the positive case for equilibrium materials, even if condensed media are employed (many typical solids have values χ (0) ∼ 10 0 -10 3 [31] and values of χ (0) ∼ 10 4 -10 5 may be observed just above the Curie point in materials that exhibit a ferroelectric phase [32,33]).On the other hand, it is essential to bear in mind that the circuit theory argument does not imply that any values in the interval −1 < χ (0) < 0 are necessarily possible for non-equilibrium materials.(Indeed, it is known, from Landau et al.'s result, that this would be an incorrect inference in the case of equilibrium materials, to which the circuit theory argument applies equally).
To determine theoretically whether values in the interval −1 < χ (0) < 0 are actually possible for non-equilibrium materials, i.e., to determine the lower permissible bound of χ (0) , it appears necessary to investigate a specific microscopic model of insulators.

C. Model
As presently conceived, materials with negative static electric susceptibility are composed of discrete entities with negative static polarizability.To begin to model condensed media consisting of such entities-be they atoms, molecules, or meta-atoms-we assume herein a set of identical, point-like, dipolarizable entities that are located at fixed positions and respond linearly and isotropically to the local electric field: in static equilibrium, the dipole moment p (j) (∈ R 3 ) of the j th entity is given by p (j) = αE (j) , where α (∈ R) is the static polarizability of each entity and E (j) (∈ R 3 ) is the ('microscopic') electric field at the point r (j) due to all sources except the j th entity itself.For a stable macroscopic medium composed of a cubic arrangement of such entities, the static electric susceptibility is also isotropic and described by scalar χ (0) (∈ R) [34].This model is, of course, precisely that which has been applied standardly, assuming α > 0, to the study of normal dielectrics with χ (0) > 0, in which context it may be traced back to the work of Mossotti [35,36], Faraday [37], and Clausius [38] (though these authors envisioned small conducting inclusions of finite, nonzero, volume rather than dipolarizable entities that were assumed point-like from the outset).The treatment of the usual case-i.e., the case α > 0 and χ (0) > 0-that we believe is most relevant to our work is that of Allen [39,40].Our task is to investigate the behavior of the model when α < 0.

D. The Drude formula and the Clausius-Mossotti equation for the case of negative static polarizability
If materials composed of entities with α < 0 are stable against self-polarization, then there is no reason to suppose that, in the weakly-interacting limit, the usual, basic expression χ (0) = Nα/ε 0 would not be as valid for the case α < 0 as it is for the case α > 0. Here, N ( 0) is the number of entities per unit volume and ε 0 is the permittivity of free space.This expression is sometimes referred to as the Drude formula for susceptibility (see, e.g., Ref. [41]) and it was assumed to be valid in the case of α < 0 by Chiao et al. [10][11][12]14]; in particular, Chiao et al. used it to arrive at their theoretical prediction of the value of χ (0) ≈ −3 × 10 −4 for pumped ammonia gas in Refs.[11,12].
The Drude formula reflects the obvious fact that, if one wishes to achieve a negative χ (0) of larger magnitude, it is natural to seek to employ entities with negative α of larger magnitude and/or arrangements of entities that are more densely packed, i.e., larger N (hence the desire for condensed media).Based on the work presented in Ref. [19], it is clear that α can be negative and the combination Nα can be of large magnitude in metamaterial systems.However, as |Nα| increases, the approximation that the entities may be treated as non-interacting-the assumption on which the derivation of the Drude formula is basedbecomes increasingly inaccurate and mutual interactions between the entities must be taken into account.Further, it must be borne in mind that, by applying the Drude formula to the α < 0 case at all, one is assuming that the χ (0) < 0 state exists; that is, one is assuming that the mutual interactions between the entities do not, in fact, make the system unstable for all values α < 0, no matter how small the value of |Nα|, in the macroscopic limit.To our knowledge, there is no definitive way to justify this assumption theoretically without taking the mutual interactions into account explicitly and confirming whether the lower limit of stability is negative definite.
Moving beyond the Drude formula, one arrives at the Clausius-Mossotti equation.The Clausius-Mossotti equation accounts for the fact that the field at a given dipolarizable entity is not, in general, just the external field, but is the sum of the external field and the field due to the polarization of all the other entities in the material.The assumptions of our model (discrete entities that are point-like, linearly-dipolarizable, at fixed positions, etc.), are exactly those often used in modern accounts of the derivation of the Clausius-Mossotti equation (see, e.g., Refs.[42,43]).The equation states that, for media with certain symmetries (including cubic crystals with one entity per primitive cell, with which we are concerned herein) [33,38,42], the static electric susceptibility is related to the static polarizability as Here, we have chosen to write the equation in terms of the dimensionless variable α = Nα/ (4πε 0 ), which will be particularly convenient later.A plot of χ (0) as a function of α according to Eq. ( 1) is included in Fig. 1.As α → 3 4π from below, the Clausius-Mossotti equation predicts that χ (0) → ∞.This is sometimes referred to as the "Lorentz 4π 3 catastrophe" and is associated with an instability towards self-polarization that has, for a long time, been considered a qualitative, if not precisely quantitative, model of the paraelectric to ferroelectric transition (see, e.g., Ref. [41]).For α > 3 4π , the Clausius-Mossotti equation formally predicts negative values of χ (0) ; however, this does not constitute a genuine prediction of a negative static electric susceptibility state.Rather, the equation ceases to be valid in this case and a more sophisticated (nonlinear) model is required to describe the self-polarized state.Accordingly, the conventional interval in which the Clausius-Mossotti equation is typically considered is 0 ≤ α < 3 4π .Since the Clausius-Mossotti equation is a useful description of dielectrics in the standard case of 0 ≤ α < 3 4π , it is natural for us to ask whether, and to what extent, it remains valid in the case of α < 0. It was noted as long ago as 1873 by Maxwell [6] that there is no step in the derivation of the Clausius-Mossotti equation that requires the static polarizability to be positive for the derivation to be valid.Therefore, if our model consists of entities that may assume negative values of α (as per the inverted media of Sanders [8,9] and Chiao et 1.The static electric susceptibility χ (0) of the material as a function of the reduced static polarizability α of the dipolarizable entities from which the material is composed, according to the model.In each of the sc, bcc, and fcc cases, the Clausius-Mossotti curve is followed, but only the system is unstable.The lower permissible bounds for bcc and fcc are identical to within the accuracy of the numerical method (orange open star), but the lower permissible bound for sc is distinct (blue open star).It is seen that the intervals of stability extend into the region of negative parameters and hence negative values of χ (0) are possible, according to the model, in each case.
The permissible lower bounds are not related to the horizontal asymptote of the Clausius-Mossotti curve (horizontal gray dotted line).
1.The limitation of the Clausius-Mossotti equation for the case of negative static polarizability Since, unlike the case of α > 0, there is no divergence of the function χ (0) for any α < 0 according to Eq. ( 1), one would expect, naively, that Eq. ( 1) should be applicable for all α < 0. If this were the case, the minimum value of χ (0) would be lim α→−∞ χ (0) ( α) = −3 (the horizontal asymptote of the Clausius-Mossotti curve as it is presented in Fig. 1).However, this is clearly at odds with the circuit theory argument, discussed previously, which has shown that any value χ (0) < −1 is unstable.
To reconcile the fact that the naive result of the Clausius-Mossotti equation for α < 0 is at odds with the circuit theory argument, one may look at the case of a sc crystal [44] of α > 0 entities to remind us of a known limitation of the Clausius-Mossotti equation that exists even within the usual interval in which it is considered to be valid, i.e., within the interval 0 ≤ α < 3 4π .The derivation of the equation would appear to apply equally to sc, bcc, and fcc crystals, and, as discussed above, Eq. ( 1) predicts that χ (0) → ∞ as α → 3 4π from below in all three cases.However, although it is not discussed explicitly in typical textbook accounts of the Clausius-Mossotti equation, it is, nevertheless, well established that a more-sophisticated analysis of the same model shows that, upon increasing α from zero, the sc case becomes unstable towards self polarization at a critical value α = α sc c + which is less than 3 4π [45] (the numerical value of α sc c + will be considered in detail below).What happens, more precisely, in the sc case is that, upon increasing α from zero, Eq. ( 1) determines correctly the associated value of χ (0) , but only until the critical value α = α sc c + < 3 4π is reached: for any value of α that is greater than α sc c + the system is already unstable and the Clausius-Mossotti curve for the sc case is truncated (as shown in Fig. 1).Accordingly, in the sc case for α > 0, the value of χ (0) is not unbounded above, but has a finite upper bound determined by evaluating χ (0) according to Eq. ( 1) at α = α sc c + .Thus, there is clear precedent for the limited applicability of the Clausius-Mossotti equation within the normally-considered interval 0 ≤ α < 3 4π , and for a finite upper bound on χ (0) , for certain crystal structures.
As is generally assumed, and we will end up essentially reaffirming in detail below, there is no such truncation of the Clausius-Mossotti curve for the bcc and fcc cases with α > 0; within the model, the upper critical values for the static polarizability are α bcc c + = α fcc c + = 3 4π and, therefore, χ (0) does indeed diverge to infinity as α → 3 4π from below.For convenience, we may refer to an instability that occurs with a divergence of χ (0) -such as that which occurs for the bcc and fcc cases on increasing α from zero-as a 'Type-I' instability, and an instability that occurs without a divergence of χ (0) -such as that which occurs for the sc case on increasing α from zero-as a 'Type-II' instability.Accordingly, to reconcile the fact that the minimum value of χ (0) = −3 that is naively predicted by the Clausius-Mossotti equation is at odds with the minimum value of χ (0) = −1 predicted by the circuit argument, we may hypothesize that, upon decreasing α from zero, there is always a Type-II instability that truncates the Clausius-Mossotti curve at or before χ (0) α = − 3 8π = −1.Sections II and III of this paper amount, essentially, to a rigorous demonstration that this is indeed what happens, and to an accurate calculation of the values of α and χ (0) at which the truncations occur, for the sc, bcc, and fcc cases.Our primary methodology is a generalization of one of the methods used by Allen-which was applied in Refs.[39,40] to the sc case with α > 0-to include also the case of α < 0.

E. The physical origin of the lower bound
We believe that, in broad terms, the physical mechanism underlying the instability is the same for the α < 0 case as for the α > 0 case (which is essentially the same for Type I and Type II instabilities).That is, the physical mechanism of the instability in the α < 0 case is essentially the same as the well-known Lorentz 4π 3 catastrophe model of the paraelectric to ferroelectric transition, which may be interpreted in the following way.It is clear that, in the absence of external sources of electric fields, the state of the system where all of the entities are unpolarized, i.e., p (j) = 0 for all j, is a state of static equilibrium of the system; for, if the polarization of each entity is exactly zero, then, in the absence of external sources, there are no electric fields that may induce polarization.The question is whether this state is a stable or unstable static equilibrium.Loosely speaking, if the entities are arranged sufficiently close together or |α| is sufficiently large, then the system may be shown to be unstable; in this case, any transient, non-zero, polarization of any of the entities would cause polarization of the other entities, which would, in turn, lead to further polarization which is, on average, of larger magnitude, etc., leading to a "runaway condition" where polarization increases, in theory, without bound.(Using a more-sophisticated, nonlinear model the self-polarization can be bounded and describe a ferroelectric or antiferroelectric state.)Conversely, if the entities are arranged sufficiently far apart or |α| is sufficiently small, then the system may be shown to be stable; in this case, any transient, non-zero polarization of any of the entities would lead to only finite oscillations of the p (j) , and these would die away to zero if any damping whatsoever were present in the system [46].

II. METHOD
We will find that any set of two or more entities (assuming the entities are located at distinct positions) has an upper critical value of α, which we denote αc + , and a lower critical value of α, which we denote αc − , such that: if αc − < α < αc + then the system is stable, and if α < αc − or α > αc + then the system is unstable.We seek to determine the values of αc + and αc − as they pertain to infinite sc, bcc, and fcc crystals, then use Eq. ( 1) to determine the associated upper and lower bounds of χ (0) .
To determine the values of αc + and αc − as they pertain to infinite crystals, our primary methodology is to consider finite crystals of a systematic range of sizes and determine the infinite limit by extrapolation.(We also summarize an alternative methodology in §III D.) This method has essentially been employed by Allen to investigate the sc case for α > 0 (see, in particular, Fig. 1 of Ref. [39]).As we will see, if the method is suitably generalized to deal also with the case α < 0, and suitable care is taken to execute and analyze the extrapolation, it can lead to definitive results with well-defined accuracy for each of the sc, bcc, and fcc cases for both α > 0 and α < 0.

A. General method for entities located at arbitrary positions
We begin by establishing the method of determining αc + and αc − for a finite number of entities at arbitrary positions, and later consider positions that represent finite sc, bcc, and fcc crystals.Consider a finite number n tot 2 of entities located at fixed position vectors r (j) (∈ R 3 ), j = 1...n tot .For now, the r (j) may be considered to be entirely arbitrary, save that no two entities are coincidentally positioned, i.e., r (j) = r (l) if j = l.As per our model, in static equilibrium the dipole moment of the j th entity, p (j) , is given by p (j) = αE (j) , where α is the static polarizability of each entity and E (j) is the electric field at the point r (j) due to all sources except the j th entity itself.In the absence of additional sources of electric fields external to the n tot dipolarizable entities themselves (an external field being superfluous to the question of the intrinsic stability of the system of entities itself), we may write that, in static equilibrium, where E (j,l) (∈ R 3 ) is the field at r (j) due to the polarization p (l) of the l th entity and l =j denotes the sum over all l = 1...n tot except for the value j.Using the standard expression for the electric field E at a displacement r from a point dipole p, viz., E = [3 (p • r) r − r 2 p] / (4πε 0 r 5 ), where r = |r|, Eq. ( 2) becomes Here, R (j,l) = r (j) − r (l) is the displacement of the j th entity with respect to the l th , and Expressing distances in multiples of a unit distance a > 0 (later, a will be the length of the sc cell edge), such that r (j) = ar ′(j) and R (j,l) = aR ′(j,l) , Eq. ( 3) may be rewritten where α = α/(4πε 0 a 3 ) is the reduced, dimensionless, static polarizability.
It is convenient to re-express the n tot vector equations of dimension three that are given by Eq. ( 4) as a single matrix equation of dimension 3n tot : or (I − αM) P = 0.
Here, I is the 3n tot × 3n tot identity matrix, 0 is the 3n tot × 1 zero column matrix, P is a 3n tot × 1 column matrix 2 , and p 3 are the Cartesian components of p (j) , and M is a 3n tot ×3n tot matrix that may be most-conveniently specified in terms of 3 × 3 sub-matrices M (i,j) , where 0 3 is the 3 × 3 zero matrix and Here, δ βγ denotes Kronecker's delta, R , and R ′(j,l) 3 are the Cartesian components of R ′(j,l) , and, in terms of the Cartesian components, R ′(j,l) = R ′(j,l) .
It is clear that the state P = 0, i.e., the state where all of the entities are unpolarized, is always a solution of Eq. ( 5) and hence always an equilibrium state of the system.This solution may be stable or unstable depending on the values of α and M, i.e., depending on the static polarizability and the relative positions of the entities.The precise criterion for stability of the system is as follows: the system is stable if all the eigenvalues of matrix (I − αM) are positive [note that the eigenvalues are real because (I − αM) is a real-symmetric matrix], and the system is unstable if any of the eigenvalues of (I − αM) are negative.
That this is the precise criterion for stability of the system may be derived using a number of different arguments.Perhaps the simplest is to note that, although we are considering entities that may exploit a source of power to maintain an unnatural state of polarization, we can, nevertheless, consider the equilibrium condition, Eq. ( 5), as arising formally from the minimization of an energy function U : R 3ntot → R whereupon Eq. ( 5) arises from the condition ∂U/∂P = 0.The matrix (I − αM) is thus identified as the Hessian matrix of the system and, according to standard results in mathematical analysis (see, e.g., Ref. [47]), the point P = 0 is a minimum of U if (I − αM) is positive definite, i.e., if all its eigenvalues are positive.If P = 0 is a minimum of the energy then it is a stable configuration of the system.Similarly, the point P = 0 is a saddle or a maximum of U if any of the eigenvalues of (I − αM) are negative.If P = 0 is a saddle or a maximum of the energy, then the system is unstable along one or more directions in P-space.
Since M is a real, symmetric, traceless matrix (of dimension 6 since we are considering n tot 2) which is not equal to the 3n tot × 3n tot zero matrix, we know that it has at least one positive (definite) eigenvalue and at least one negative (definite) eigenvalue, and the condition for stability may be restated most simply as follows: the system is stable for values of α such that αc − < α < αc + , where αc − = 1/λ min < 0 and αc + = 1/λ max > 0. Here, λ min and λ max are the minimum and maximum eigenvalues of M respectively.Conversely, if α < αc − or α > αc + then the system is unstable.
Thus, the question of stability of a given system becomes essentially a question of constructing the matrix M for that system (which depends only on the relative positions of the entities) and calculating its maximum and minimum eigenvalues.This stability criterion is a generalization-to include also the α < 0 case-of that presented, and derived by a somewhat similar argument, for the α > 0 case in Refs.[39,40].

B. Application to arrays forming finite sc, bcc, and fcc crystals
Having established the stability criterion for entities located at arbitrary positions, we now proceed to specify explicitly the locations of the entities for sc, bcc, and fcc finite crystals used in the study.At this stage, it is convenient to switch from labeling the entities with a single index j (with position vectors r (j) , etc.) to labeling them with an ordered triple of integers (u 1 , u 2 , u 3 ).In this notation, finite crystals may be created via the usual approach [48] of locating entities at the sets of (reduced) position vectors ), and ãfcc 3 = (1/2 1/3 ) (x 1 + x2 ).Each set S n describes a finite crystal consisting of an n × n × n array of entities (total number of entities n tot = n 3 ) with the overall shape of a rhombohedron (more specifically, a cube in the sc case).
We have chosen normalization factors, 1/2 2/3 and 1/2 1/3 for the bcc and fcc cases respectively, such that the primitive cells (and not the conventional cells) are of unit volume.This means that the number density is the same in each of the sc, bcc, and fcc cases and given by N = 1/a 3 , which is convenient for our purposes.In particular, the expression for α used in this section, α = α/(4πε 0 a 3 ), is consistent with the expression α = Nα/(4πε 0 ) used in Section I C and Eq. ( 1), and the resulting values of αc ± for the sc, bcc, and fcc cases may be inserted into Eq.( 1) and compared in a like-for-like fashion.
We used Python to generate the M-matrices, calculate λ min and λ max , and hence determine the critical values of the static polarizability αc + , n and αc − , n for a given n × n × n crystal described by a given set S n .The Python code is included in the Supplemental Material [49].
The value n max = 27 was used as it was found to be the largest value for which our code could be run on the system with 64 GB RAM that we had readily available and, as shown below, it already provides results that are sufficiently conclusive for our purposes.
We may note that the method takes into account all interactions between every pair of entities in the crystal (and does not, for example, assume only nearest neighbor interactions).
To determine the critical values of the static polarizability in the macroscopic limit, i.e., in the limit of infinite sc, bcc, and fcc crystals we calculate the limiting values αc + , ∞ = lim n→∞ αc + , n and αc − , ∞ = lim n→∞ αc − , n by extrapolation, as detailed below.

A. Critical values of the reduced static polarizability and static electric susceptibility
Raw data for the values of αc − , n and αc + , n in each of the sc, bcc, and fcc cases are given in Table S1 of the Supplemental Material [49].Plots of αc − , n and αc + , n for n = 2, 3, ..., n max are shown in Fig. 2. It is seen from the plots that: (1) In each of the sc, bcc, and fcc cases, the values of αc + , n decrease with increasing n and the values of αc − , n increase with increasing n.
(For the case of α sc c + , n , the decreasing behavior was essentially noted in [39], although there it was characterized as an increase in the associated eigenvalue, which amounts to the same thing, since αc + = 1/λ max .)Thus, the interval of stability of α is reduced in both its positive and negative extents as n increases.This makes intuitive physical sense; it would seem reasonable to expect that the addition of more dipolarizable entities can, loosely speaking, only serve to increase the overall amount of mutual interaction within the system, causing a greater tendency towards instability.( 2) In all cases, the sequences αc ± , n , n = 2, 3, 4, ...  shown in Fig. 3.The plots indicate also how the limiting values αc ± , ∞ may be determined by extrapolation.It can be seen that the plots are approximately linear in each case; this enables the best estimates for the values of αc ± , ∞ = lim n→∞ αc ± , n to be determined as the vertical intercepts of linear fit lines.Since the data becomes linear to an increasing level of accuracy as n is increased, we use the two largest-n data points available, n = n max − 1 = 26 and n = n max = 27, to construct the extrapolation line and determine what we believe to be our best estimates of the limiting values using this method.Explicitly, the formula for the extrapolated values is thus The values obtained by this method are listed in Table I, where they are quoted to three significant figures.
As discussed above, we expect the Clausius-Mossotti equation to remain valid for values of α that are within the stable interval αc − < α < αc + .Therefore, the critical values of the static electric susceptibility χ c ± may be determined by inserting the values of αc ± , ∞ into Eq.( 1).The values of χ (0) c ± thus obtained are listed also in Table I to three significant figures.The results of the study are summarized in Fig. 1.

B. Accuracy of the numerically-determined critical values
There are a number of ways by which we may consider the accuracy of the results determined via the above method.
1. Accuracy via comparison against assumed exact values for the cases of bcc and fcc with α > 0 It appears reasonable to assume that the true model values (by which we mean the values predicted by the model, if we were able to solve it exactly as opposed to numerically) of α bcc c + , ∞ and α fcc c + , ∞ are exactly 3 4π , as per the vertical asymptote in the Clausius-Mossotti expression, Eq. (1); our numerical results for these values are α bcc c + , ∞ = 0.238752 and α fcc c + , ∞ = 0.238756 to six decimal places, which approximate 3 4π to an accuracy of 0.008 % and 0.010 % respectively.Based on this, it appears reasonable to argue that, in general, the uncertainty in determining the values of αc ± , ∞ via our method is probably 0.01 % to the nearest order of magnitude (when using n max = 27); that is, we may assume an uncertainty of the order 0.01 % also on our values of α sc c + , ∞ , α sc c − , ∞ , α bcc c − , ∞ , and α fcc c − , ∞ .Using the standard method for propagation of uncertainties, uncertainties of the order 0.01 % on α sc c + , ∞ , α sc c − , ∞ , α bcc c − , ∞ , and α fcc c − , ∞ , lead, via Eq.( 1), to uncertainties of the order of 0.01 % also on the values of χ ).Therefore, we may conclude that the values listed in Table I are almost certainly accurate as-quoted to three significant figures, and, if desired, values quoted to a larger number of significant figures with uncertainties of the order of 0.01 % may be readily generated from the raw data (provided in the Supplemental Material [49]).In particular: to a larger number of significant figures, we find χ  1), be exactly − 3 8π .Our numerical results for these values are α bcc c − , ∞ = −0.119382and α fcc c − , ∞ = −0.119381 to six decimal places, which may be seen to approximate − 3 8π to an accuracy of 0.013 % and 0.012 % respectively.This appears consistent with our previous assumption that the uncertainty on all the numerical values of αc ± , ∞ is of the order 0.01 %.

Accuracy via inductive reasoning
An alternative way to assess the accuracy of our numerical results is as follows.Noting, from Fig. 2 and the raw data tabulated in the Supplemental Material [49], that the sequences αc + , n , n = 2, 3, 4, ..., 27, decrease monotonically with increasing n in each of the sc, bcc, and fcc cases, and assuming that this trend continues for all n (i.e., applying 'inductive reasoning'), then any particular value of αc + , n for a given crystal structure provides an upper bound on αc + , ∞ for that crystal structure (note that we are referring here to an upper bound on the numerically-determined value of αc + , ∞ , which is itself an upper bound on the stable value of α).Accordingly, the most stringent upper bounds on the values of αc + , ∞ that may be determined from the data in this way are given by αub c + , ∞ = αc + , nmax (with, in our case, n max = 27).
To determine lower bounds on the numerically-determined values of αc + , ∞ , we first recall that our extrapolation to determine the best-estimate values of αc ± , ∞ , as per Fig. 3 and Eq. ( 9), used a straight line fit through the points (1/n 2 , αc + , n ).Here, the exponent of two in 1/n 2 was chosen 'by hand' to produce the most-linear plot.If, instead of an exponent of two, we choose an exponent of one and, hence, plot instead the points (1/n, αc + , n ), we may observe that the graph is increasing and convex in each of the sc, bcc, and fcc cases (see Fig. S1 of the Supplemental Material [49]).The choice of an exponent of one is somewhat arbitrary; the requirement is simply that a convex graph is produced.Again, assuming this trend remains true, i.e., the plot remains increasing and convex, not just for n = 2...27 but for all n, then the vertical intercept of a straight line through any two data points provides a lower bound on αc + , ∞ .Accordingly, the most stringent lower bounds on the values of αc + , ∞ that may be determined in this way from the data available are given by the vertical intercepts of the straight lines through the points (1/n, αc + , n ) with n = n max − 1 = 26 and n = n max = 27.Explicitly, the formula for this procedure is thus Similarly, with regard to αc − , ∞ , we may observe that the sequence αc − , n , n = 2, 3, 4, ..., 27 increases monotonically for each type of crystal structure and we may argue that lower bounds on the values of αc − , ∞ are given by αlb c − , ∞ = αc − , nmax in each case.To determine upper bounds on the values of αc − , ∞ , we observe that plots of the points (1/n, αc − , n ) are decreasing and concave in each case, with the caveat that only points for which n 5 are included in the sc case (see Fig. S1 of the Supplemental Material [49]).Therefore, we may argue, analogously to above, that the most stringent upper bounds on the values of αc − , ∞ that may be determined by this method from the data available are given by the vertical intercepts of the straight lines through the points (1/n, αc − , n ) with n = n max − 1 = 26 and n = n max = 27.Explicitly, the formula for this procedure is thus In this way we find, to four significant figures: α sc c + , ∞ = 0.1868 (More specifically, slab-like square cuboids in the sc case.) • 'Parallelepiped needles', formed from entities at the sets of position vectors (More specifically, needle-like square cuboids in the sc case.) • Spheres, formed from entities at the sets of position vectors In each case, if the values of αc + , n and αc − , n for set S n are calculated for a range of n, the limits αc ± , ∞ = lim n→∞ αc ± , n are determined by extrapolation, and the accuracy is determined by one or other of the methods above-i.e., if the same process is carried out for these alternatively-shaped crystals as was carried out for the rhombohedral crystals abovethen the same values of α sc c ± , ∞ , α fcc c ± , ∞ , and α bcc c ± , ∞ , and hence the same values of χ are found, within the accuracy of the method (i.e., the same results as already reported in Table I).
Of course, this approach does not rule out the possibility that some other overall shape of crystals that we have not checked explicitly-say, ellipsoidal-may somehow give different values, but we nevertheless consider it appropriate to reasonably conclude that the method appears to provide results that are independent of the sample shape.
D. An alternative methodology (yields the same results) We may refer to the method employed hitherto in this paper as the 'finite crystal method' since the macroscopic limit is considered by extrapolating the results for finite crystals of increasing size.An alternative methodology involves assuming an infinite crystal from the outset, which we may refer to as the 'infinite crystal method'.Both methods were considered and applied by Allen to investigate the sc case with positive polarizability in Ref. [39].
Each method has certain advantages and disadvantages, but we believe that, overall, the finite crystal method provides the most direct and rigorous route to determine the values of αc ± with well-defined accuracy in the macroscopic limit, and hence it is the method we reported in detail above.However, we have also carried out a detailed analysis of the infinite crystal method as applied to sc, bcc, and fcc crystals with emphasis on the negative static polarizability case [50]; in all cases, we find excellent agreement with the above-stated results, providing a useful cross-check and validation.
One favorable feature of the infinite crystal method is that, once the location in reciprocal space of the extremal eigenvalues has been established for a given crystal structure, the values of αc ± (which are, necessarily, the macroscopic values) may be expressed as infinite lattice sums.For example, in the sc case we may write and where u ′ denotes the sum over all triples of integers u = (u 1 , u 2 , u 3 ) except (0, 0, 0) (corresponding to the interaction of a given dipole with all other dipoles in the infinite crystal).
We are not aware of any closed form expressions for these two particular sums, but they are absolutely convergent and may be readily approximated either via brute force summation with u 1 , u 2 , u 3 = −N sum ...N sum for some large N sum , by evaluating the sum for a range of N sum and extrapolating N sum → ∞, or by using the Ewald summation method.To three significant figures, the values of the sums are α sc c + = 0.187 and α sc c − = −0.103,which agree with the results found previously for the finite crystal method recorded in Table I.
The reciprocal λ sc max = 1/ α sc c + = 5.35 was considered explicitly by Allen and is identified as the maximum point of the graph in Fig. 3 of Ref. [39].The lower critical values were not of interest to Allen, but we may see that the reciprocal λ sc min = 1/ α sc c − = −9.69 is the minimum point of the graph in Fig. 3 of Ref. [39].(We are not aware of any general theorem that states the extremal eigenvalues must lie along the edges of the irreducible Brillouin zone, as Allen appears to assume, but our work [50] indicates that this is indeed true for the sc, bcc, and fcc cases we have considered.) The lattice sums of Eqs. ( 12) and ( 13) appear also in other physical systems and have, long ago, been evaluated in other contexts.For example, they appear in the related problem Our initial motivation for carrying out the work presented herein was a concern that mutual interactions between entities with α < 0 might, collectively, wipe out the possibility of χ (0) < 0 materials in the macroscopic limit: that is, we were concerned that, in the limit of large particle numbers, the lower critical value of α would converge to zero and, accordingly, the lower permissible bound of χ (0) would be zero in all cases.(The experiments reported in Ref. [19] were not sophisticated enough to rule out this possibility.)Perhaps most importantly, then, we have seen that this does not happen; the lower permissible bounds of χ (0) are negative definite in each case.Therefore, we conclude that negative static electric susceptibility values are indeed possible in non-equilibrium crystals, according to the model.(We note that this conclusion is very different from that of the circuit theory argument, which also showed that values χ (0) < −1 are impossible whether the material is in equilibrium or not but did not imply that any values in the interval −1 < χ (0) < 0 are necessarily possible.) We see that the lower permissible bound depends on the crystal structure: it is different for sc than for bcc and fcc.Thus, the question of the lower permissible bound for a given crystal structure appears to be, in general, a highly non-trivial question which can probably only be addressed definitively by the type of methods considered herein.In particular, we can assert that there is no structure-independent 'short cut' derivation of the lower permissible bound.
Although the value χ = −1.00,and the assumption of point-like dipolarizable entities means that the results of the study may be inaccurate in relation to any real system, our results suggest that, if one wishes to create an isotropic material with the lowest possible value of χ (0) , a bcc or fcc structure, as opposed to a sc structure, is likely to be preferable.
It is clear that the case of χ (0) < 0 is very different from the usual case of χ (0) > 0.
Whereas the permissible value of χ (0) may be unbounded above for certain crystal structures (e.g., bcc and fcc with one entity per primitive cell), it is always bounded below by the value −1 according to the circuit theory argument, and, as we have seen, may be even more limited in its negative extent in certain crystal structures (e.g., sc with one entity per primitive cell).
Nevertheless, in finding that the permissible interval may extend all the way down to the circuit theory limit of −1 for certain condensed media, we have confirmed that condensed media are, indeed, capable of exhibiting χ (0) < 0 values with magnitudes that are ∼ 10 3 times greater than those proposed initially in gaseous systems, even when mutual interactions between the dipolarizable entities are taken into account.We believe that this increase is likely to be sufficient to enable the remarkable potentialities of χ (0) < 0 materials-such as new forms of charged particle traps-to become practically feasible. FIG.

FIG. 2 .
FIG.2.The positive and negative critical values of the reduced static polarizability, αc + , n and αc − , n respectively, for sc, bcc, and fcc finite crystals as a function of n, for n = 2, 3, 4, ..., 27 (total number of entities n tot = 2 3 , 3 3 , 4 3 , ..., 273 ).The vertical lines denote the intervals of reduced static polarizability αc − , n < α < αc + , n that are stable in each case.For clarity, the data sets for bcc and fcc have been displayed with artificial offsets in the horizontal direction of n = 0.2 and n = 0.4 respectively.
appear to converge to finite, non-zero values.(We already know, by the very fact that the dielectric state exists for the case of positive static polarizability within this model, that the sequences αc + , n , n = 2, 3, 4, ... do not converge to zero for large n, and the qualitative behavior of the magnitude of sequences αc − , n , n = 2, 3, 4, ... appears to be very similar in this respect.)(3) In both the positive and negative cases, it appears that the values of α bcc c, n and α fcc c, n become increasingly coincident with each other as n increases, whereas the values of α sc c, n remain significantly different.These aspects of the data are illuminated further by plots of the points (1/n 2 , αc ± , n ), as

TABLE I .
Extrapolated critical values of the reduced static polarizability αc ± , ∞ , and the associated critical values of the static electric susceptibility χ ± , for sc, bcc, and fcc infinite crystals, stated to three significant figures.Where available, the corresponding exact values (expected or inferred)