L F Pašteka

Trygve Helgaker
email: t.u.helgaker@kjemi.uio.no

Trond Saue
email: trond.saue@irsamc.ups-tlse.fr

Dage Sundholm
email: sundholm@chem.helsinki.fi

Hans-Joachim Werner
email: werner@theochem.uni-stuttgart.de

M Hasanbulli
email: mustafa.hasanbulli@gmail.com

J Major

Peter Schwerdtfeger
email: p.a.schwerdtfeger@massey.ac.nz

H.-J Werner

Atoms and molecules in soft confinement potentials

published or not. The documents may come

I. INTRODUCTION

The confined hydrogen atom was introduced more than 80 years ago by Michels et al. [1] in 1937, who studied a system consisting of one hydrogen atom at the center of an impenetrable spherical confinement potential. Their goal was to observe how atomic properties evolve as a function of compression -that is, with changing the confinement radius r 0 . Since then, many authors have investigated similar model systems consisting of atoms or molecules confined by impenetrable or partially penetrable walls of different geometrical shapes such as spherical, paraboloidal, or prolate spheroidal walls [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. An example of such a system is the non-relativistic, artificially bounded harmonic oscillator, enclosed between potential walls. This model has been successfully applied to problems such as the fundamental mass-radius relation for white dwarf stars [START_REF] Auluck | Proc. Nat. Inst. Sci. India[END_REF], the rate of escape of stars from galactic and globular clusters [START_REF] Chandrasekhar | [END_REF], the role of the symmetrically bounded linear harmonic oscillator in the theory of the specific heat of solids [19], second-order phase transitions [20], energy levels and transition probabilities for a bounded linear oscillator [21], anharmonic effects in solids [22], magnetic properties of metallic solids [23], and nuclear shell models [24]. Similar model systems have been employed in various research fields where the effects of pressure on energy levels [1,8,12] and properties such as polarizabilities [1,9], hyperfine splittings [7][8][9]11], nuclear magnetic shielding constants [9], hyperfine interaction energies [10] and electron (de)localization [25] have been of interest. We note that quark confinement is at the core of proton stability [26], a necessary condition for matter as we know it to exist. More detailed accounts of the work on confined systems are found in Refs. [27][28][29][30][31].

Over the years, a wide range of analytical and numerical methods have been employed to obtain solutions for confined systems. Many of the early studies modified and adapted the well-known analytical solutions for the free hydrogen atom -see, for example, Refs. [2,32,33]. Already in 1938, Sommerfeld and Welker [2] presented the formal solution to the problem proposed by Michels et al. [1] in terms of confluent hypergeometric functions. Since the lack of computational resources made it impossible to obtain accurate energy eigenvalues, they were calculated using analytical expansions of the hypergeometric functions. Sommerfeld and Welker also identified the critical cage radius r HW c at which the energy of the lowest eigenvalue becomes zero. Below this radius, the electron of the hydrogen atom is bound by the confinement potential rather than by the Coulomb potential. In 1946, de Groot and ten Seldam [3] studied energy eigenvalues of excited states, considering the problem with non-zero angular momentum. In addition, these authors solved the problem for total energy E > 0 and were the first to realize that, for small values of r 0 , the kinetic energy of the system is higher than the Coulomb potential energy.

After Sommerfeld and Welker's [2] reformulation of the problem, many authors revisited this model. Suryanarayana and Weil [7] formulated the wave functions in terms of the confluent hypergeometric functions and calculated energy eigenvalues numerically to a few significant figures. Both Goldman and Joslin [34] and Chuu et al. [35] used a modified version of Kummer's differential equation and represented the formal solutions in terms of Whittaker functions, computing energy eigenvalues from a truncated series expansion of the Whittaker functions with many terms. In 2005, Burrows and Cohen [36,37] investigated the model using a combination of group theory and algebraic methods. Recently, in 2007, Aquino et al. [38] showed that it is possible to obtain the energy for the confined hydrogen atom to a very high accuracy, with up to hundred decimal digits. In addition to energy eigenvalues for the ground state and many excited states, they computed expectation values r -1 , r and r 2 , hyperfine splitting and magnetic screening constants, polarizabilities in the Kirkwood approximation, and pressure as a function of the confinement radius.

Interest in the confined hydrogen atom [39] also kindled the interest in confined multi-electron systems. The helium atom is the simplest system where electron correlation can be investigated. The first results on the confined helium atom originated in 1952 by ten Seldam and de Groot [40]. Based on a trial wave function proposed by Hylleraas in 1929 [41], they performed a variational calculation on the confined helium system, calculating the pressure by differentiating the energy with respect to the volume. In a second publication [42], they calculated the polarizability of helium and observed the expected decreasing trend with respect to decreasing confinement radius.

Since ten Seldam and de Groot's attempts, a vast amount of articles have focused on energy shifts of ground and excited states and various electronic-structure properties of the confined helium atom. The most up-todate results were reported in 2010 by Montgomery Jr. et al. [43], who solved the time-independent Schrödinger equation using first-order Rayleigh-Schrödinger perturbation theory and then extended their calculations using fifth-order variational perturbation theory. Using density-functional theory, Waugh et al. [44] looked into the variation of the static dipole polarizability and hyperpolarizability with the strength of confinement for the helium atom. Most recently, de Oliveira Batael and Drigo Filho investigated H 2 in an impenetrable prolate spheroidal box [45].

Most studies presented so far have considered Coulomb plus hard-wall potentials, which are ill suited for traditional quantum-chemistry methods using basis functions such as Gaussians, which extend to infinity. In this article, we focus on a soft spherical confining potential of the form (r/r 0 ) N with stiffness parameter N ∈ N and confinement radius r 0 . We apply such a confinement potential to the non-relativistic hydrogen, helium, carbon, and potassium atoms and to the CH 4 molecule in the center of the spherical confinement, noting that the relativistic Coulomb plus confinement potentials entails additional difficulties [46]. For confined hydrogen, a more general potential including this form was introduced by Diercksen and co-workers [47,48] and more recently investigated by Katriel and Montgomery, who also discussed the virial theorem and convergence towards the hardwall limit [49]. The harmonic-oscillator case (N = 2) has been used extensively for Coulombic systems (see Refs. [48,50,51] and references therein), while the linear case (N = 1) has been used in physics for states of quarkonia [52]. Interesting applications are, for example, the stability of doubly negatively charged ions in external confinement potentials [53]. We should also mention a very recent study of Rahm et al. on isotropic compression of single atoms in a non-reactive neon-like environment using a pressure polarizable continuum model [54].

Like similar confining potentials, the soft confinement potential (r/r 0 ) N allows us to study the effects of the spatial restrictions on a system by changing the confinement radius r 0 . In addition, by varying the stiffness parameter N , we can observe the transition from a system placed in a 'soft' spherical box to one that is confined in an hardwall spherical box. Importantly, unlike hard-wall confined many-electron systems, a soft-wall potential does not possess a discontinuity [49], leading to eigenfunctions defined over the whole complex space. This lack of discontinuity simplifies its use with existing numerical algorithms and with standard Gaussian-based electronicstructure software packages enormously, opening the way for quantum-theoretical studies of atoms or molecules in relatively stiff potentials, with N 0. Many-electron systems in a soft confinement potential have not been investigated in great detail except for the harmonic case, which is the reason for the present study.

II. SOFT-WALL CONFINEMENT

A. Confinement potentials

The Hamiltonian of a spherically confined n-electron atom with nuclear charge Z is given by (in atomic units)

H c = n i=1 - 1 2 ∇ 2 i - Z r i + V conf (r i ) + n i=1 n j>i 1 r ij (1)
where r i = |r i | and r ij = |r i -r j |, and where V conf (r i ) represents the spherical confining potential. The archetypal choice for V conf is a spherical box of confinement radius r 0 with an impenetrable wall (the hard-wall potential) [1]:

V HW (r) = 0 if 0 ≤ r < r 0 , +∞ if r ≥ r 0 . (2)
Because of the discontinuity of the potential at r 0 , any eigenfunction ψ(r 1 , . . . r n ) of the Hamiltonian Eq. ( 1) with V conf (r i ) = V HW (r i ) must vanish when r i = r 0 for some coordinate i, satisfying the Dirichlet boundary conditions ψ(r 1 , . . . r n ) = 0 at the boundary surface.

Here, we consider instead the soft-wall (or polynomial) confining potential, introduced early on by Diercksen and co-workers [47,48,55]. For a multi-electron system, it is given by

V N (r i ) = n i=1 r i r 0 N , (3)
where r 0 is the radius of the spherical box and N ∈ N the stiffness parameter (the nuclei are confined by placing the center of mass at the origin of the sphere). By increasing N , the potential becomes steeper at r 0 , thereby making the wall stiffer and less penetrable. Considering the limit N → +∞, it is readily seen that V N (r) → V HW (r) point-wise as shown in Figure 1. Yet, in contrast to V HW (r), the potential V N (r) does not possess a discontinuity at r 0 , thereby removing the necessity to employ the Dirichlet boundary condition. Still, it can be expected that any eigenfunctions of the Hamiltonian Eq. ( 1) with V conf (r i ) = V N (r i ) will decay quite rapidly for |r i | > r 0 .

B. Virial theorem for confined atoms

The virial theorem for homogeneous potentials was considered by Fock already in 1930 [56]. Following his scaling argument, we here consider a many-electron system described by the n-electron Hamiltonian of Eq. ( 1) written in the form

H = T + V C + V N , ( 4
)
where T is the kinetic-energy operator T = -

1 2 i ∇ 2 i , V C is
V C = -i Zr -1 i + i>j r -1 ij = V 1 +V 2
, and V N is the soft confinement potential defined in Eq. ( 3). Note that the Coulomb and confinement potentials are homogeneous Laurent polynomials of degrees -1 and N , respectively.

The n-electron ground-state energy

E = Ψ|T + V C + V N |Ψ = T + V C + V N (5)
is stationary with respect to all norm-conserving variations in the exact n-electron ground-state wave function Ψ. Consider, in particular, the norm-conserving uniform scaling of all coordinates:

Ψ α (r i ) = α 3n/2 Ψ(αr i ) (6)
such that Ψ 1 = Ψ. Some elementary considerations show that the expectation value for the scaled wave function is given by

Ψ α |H|Ψ α = Ψ|α 2 T + αV C + α -N V N |Ψ . (7)
From the Hellmann-Feynman theorem, we then obtain

d Ψ α |H|Ψ α dα α=1 = Ψ|2T + V C -N V N |Ψ . (8)
Using the stationarity of the energy for α = 1 to set the result equal to zero, we arrive at the virial relation for soft confinement:

2 T + V C = N V N . (9)
Using scaling arguments, Katriel and Montgomery established this expression for one-electron systems [49]. In the absence of a confining potential (i.e., when N = 0), it reduces to the standard virial relation 2 T + V C = 0 for atoms. Differentiating the ground-state energy with respect to the confinement radius r 0 and invoking the Hellmann-Feynman theorem again, we obtain

dE dr 0 = Ψ ∂V N ∂r 0 Ψ = - N r 0 V N . (10)
For soft confinement, we may therefore write the virial relation in the alternative form

2 T + V C = -r 0 dE dr 0 . (11)
Since the stiffness parameter N does not appear explicitly in this expression, it must hold also for hard confinement, in the limit N → ∞. Indeed, the virial theorem in this form was established for hard-wall confinement by Fernández and Castro [57].

Using the virial relation in Eq. ( 9) to eliminate V N from the expression for the total energy in Eq. ( 5), we obtain

E = N + 2 N T + N + 1 N V C . (12)
Letting N → ∞ and comparing with the energy expression in Eq. ( 5), we find that the confinement potential does not contribute directly to the energy in the hardwall limit lim

N →∞ V N = 0 , (13)
only indirectly by setting up the boundary conditions, in agreement with the energy expression E = T + V C for hard-wall confinement.

Let us now consider the relation between the Coulomb energy and the kinetic energy for small confinement radii. Dividing the virial relation in Eq. ( 9) by V C , we obtain the equivalent relation

2 T V C + 1 = N V N V C . ( 14
)
Letting r 0 → 0 on both sides for fixed N , we find that if

lim r0→0 T V C = ∞ ⇐⇒ lim r0→0 V N V C = ∞ . (15)
The kinetic energy thus dominates over the Coulomb energy if (and only if) the confinement potential dominates over the Coulomb energy. This result holds for soft confinement potentials with N > 1. Letting N → ∞, we find that it holds also for the hard-wall confinement. This behavior of the kinetic and Coulomb energies for small confinement radii is also suggested by the behavior of the expectation values of T and V C upon scaling of the ground-state wave function Ψ α . From the expression in Eq. (7), we find that the kinetic energy increases quadratically and the Coulomb energy only linearly as the wave function Ψ α contracts uniformly with increasing α. Since the ground-state wave function may contract in a non-uniform manner as the confinement radius becomes small, the ground-state kinetic and Coulomb energies may behave differently.

III. THE SOFT-WALL CONFINED HYDROGEN-LIKE ATOM

A. Theory

For a confined hydrogen-like atom with nuclear charge Z, we must solve the radial Schrödinger equation

- 1 2 d 2 dr 2 + ( + 1) 2r 2 - Z r + V conf (r) P n (r) = E n P n (r).
(16) with boundary conditions P n (r) → 0 for r → 0 and r → ∞. In the absence of confinement, solutions may be expressed in terms of confluent hypergeometric functions as

P n (r) ∼ x +1 F ( + 1 -n; 2 + 2; x)e -x/2 , (17)
with x = 2kr, where k 2 = -2E n and the quantum number n = Z/k ∈ N is introduced to have bound solutions. Under confinement, it is convenient to keep n, but it now serves as a counting number [START_REF] Scherbinin | Physics of Clusters[END_REF][START_REF] Ferreyra | [END_REF], as will be further elaborated in the following. No analytical solution with the potential V conf (r) = V N (r) over the whole domain is known, and the well-known n 2 orbital degeneracy of the E n levels for a given principal quantum number n and different < n values present for V conf (r) = 0 is lifted under confinement [3]. More specifically, for a given principal quantum number n > 0 of the hard-wall confined hydrogen-like atom (V conf (r) = V HW (r)), states of high become stabilized over states of low upon confinement, leading to several excited-state level crossings (accidental degeneracies) at small confinement radii r 0 [39,[START_REF] Scherbinin | Physics of Clusters[END_REF][START_REF] Ferreyra | [END_REF][60][61].

For the soft-wall confined hydrogen atom (V conf (r) = V N (r)), the behavior is similar, as shown by treating V N (r) as a perturbation to the free hydrogen atom at large r 0 . For fixed n, we must then consider the expectation values r N nl for different values of < n. We note that since P n (r) vanishes as r +1 as r → 0 the expectation values r M nl are divergent for real M ≤ -2 -3, as illustrated in Figure 2. A closed formula for mean values of all powers in r in hydrogenlike atoms has been given recently in a seminal paper by Suslov and Trey [62], which also contains a rich bibliography of work in this area (see also [63]). Here we shall rather use the recurrence relations developed independently by Pasternack [START_REF] Pasternack | Proc. National Acad. Sci[END_REF] and Kramers [START_REF] Kramers | Die Grundlagen der Quantentheorie/Quantentheorie des Elektrons under der Strahlung[END_REF]. Whereas Pasternack exploited properties of hypergeometric functions, Kramers multiplied the radial equation, Eq. ( 16), (with V conf (r) = 0) by

r M +1 dP n dr - 1 2 (M + 1) r M P n , (18)
and obtained, after radial integration (with suitable integration by parts) and inserting the energy expression for bound states, the Kramers-Pasternack recurrence re- lation (in atomic units)

(M + 1) Z 2 n 2 r M n -(2M + 1)Z r M -1 n + 1 4 M (2 + 1) 2 -M 2 r M -2 n = 0 , (19)
(Z > 0, ≥ 0, n ≥ 1, r M n,l > 0) with the following initial terms for upward and downward recursion in M , respectively

r 0 = 1, r n = 1 2Z 3n 2 -( + 1) , (20)
r -1 n = Z n 2 , r -2 n = 2Z 2 n 3 (2 + 1)
, and the identity

r M -1 = (2 + M + 1)! (2 -M )! n 2Z 2M +1 r -M -2 . ( 21
)
for 0 ≤ M ≤ 2 . The expectation value r -1 follows directly from the virial theorem, and the r -2 follows immediately from (21). From these relations, the strict inequality

r N n, -1 > r N n ( 22
)
for a fixed Z value follows. Surprisingly, a proof for this inequality is not available in the literature. To fill this gap, we provide a proof by induction in Appendix A. As a result, for the confining potential Eq. ( 3) with a finite radius r 0 1 and stiffness parameter N ∈ N, we obtain the following strict inequality to first order in perturbation theory,

E n, -1 (r 0 , N ) > E n (r 0 , N ), (23)
which holds in the hard-wall limit as well [START_REF] Scherbinin | Physics of Clusters[END_REF].

We also introduce a second inequality

Z r N +1 n -n 2 r N n > 0 ∀n, < n, N ∈ N 0 (24) n 2 r -N -1 n -Z r -N n > 0 ∀n, < n, N ∈ N (25)
and for the latter expression we are limited to N < 2 + 2 (for N =0 the latter equation becomes exactly zero). Note that the r-expectation values are Z-dependent. Further, if (24) holds, (25) follows from the identity (21). A proof of eq.( 24) is provided in Appendix B. In first-order perturbation theory for finite r 0 and N ∈ N, we therefore get

E n (r 0 , N + 1) > n 2 Zr 0 E n (r 0 , N ). ( 26
)
With increasing stiffness, this implies an increasing effect on the energy levels with higher principal quantum number and smaller confinement radius, as is intuitively clear. The relativistic formulation of the Kramers-Pasternack recurrence relation is far more complicated [START_REF] Suslov | [END_REF]67] and similar inequalities still need to be explored.

It is interesting to note that the stabilization of states of high angular momentum over those of low angular momentum is observed also for atoms and molecules in magnetic fields [68]. In a non-relativistic treatment, the magnetic field affects the energy levels in two ways: paramagnetically through the spin and orbital Zeeman operators (splitting energy levels, favoring the states with the lowest component of angular momentum in the field direction) and the diamagnetic operator (raising the energy of all states). The diamagnetic operator corresponds to a harmonic confinement potential in the plane perpendicular to the field direction, raising the energy of all states but least for the radially more compact highangular-momentum states.

Before closing this section, we note that if we follow the approach of Kramers [START_REF] Kramers | Die Grundlagen der Quantentheorie/Quantentheorie des Elektrons under der Strahlung[END_REF], starting from the radial equation, Eq. ( 16), with a soft confinement potential

V conf (r) = V N (r)), we obtain -2E(M + 1) r M -(2M + 1)Z r M -1 (27)
+ M 4 (2 + 1) 2 -M 2 r M -2 + (2M + 2 + N ) r N 0 r M +N = 0 .
Again, setting M = 0 we obtain the energy expression

E = 1 2 V C + N + 2 2 V N (28)
which can also be obtained from the virial relation, Eq. ( 9).

B. Z Scaling

Here we analyze the Z scaling for hydrogen-like atoms in a soft confinement potential. Hylleraas used coordi-nate scaling r → r/Z to show that hydrogenic eigenvalues and eigenstates for Z > 1 have simple relations to the Z = 1 case [69,70]:

E Z = Z 2 E Z=1 ; ϕ Z (r) = Z 3/2 ϕ Z=1 (Zr) . (29)
Introducing this scaling for the confined hydrogen atom in Eq. ( 16), we obtain

Z 2 - 1 2 d 2 dr 2 + ( + 1) 2r 2 - Z r + Z -N r r 0 N P n (r/Z) = E n P n (r/Z) . (30)
To connect one-electron atoms of different nuclear charge, we must also scale the confinement radius r 0 :

r N 0 → r N 0 /Z N +2 =⇒ r 0 → r 0 /Z 1+ 2 N , (31)
which in the hard-wall limit (N → ∞) gives r 0 → r 0 /Z. The energy of a hard-wall confined one-electron atom of charge Z and confinement radius r 0 /Z is therefore given by Z 2 E n, where E n, is the energy of a one-electron atom of unit nuclear charge and confinement radius r 0 , as already pointed out by Patil and Varshni [71]. Hence, without loss of generality, we can restrict our study of one-electron atoms to the hydrogen atom, with Z = 1. From the recursion formula Eq. ( 19), we immediately see that r N n ;Z = Z -N r N n ;Z=1 , ∀n ∈ Z, < n. These observations further allow the extension of the 1/Z expansion for the ground-state energy of the free helium isoelectronic sequence to the confined case, as will be discussed in Section IV B.

C. Asymptotic behavior

The asymptotic behavior of the radial wave functions P n (r) for confined hydrogen-like atoms as r → 0 is identical to that of the free atom, that is, P n (r) vanishes as r +1 . On the other hand, the asymptotic behavior of the radial functions P n (r) as r → +∞ is expected to be markedly different and needs to be analyzed for future basis-set studies. In the asymptotic limit and with N > 0, the term in brackets of Eq. ( 16) is dominated by the kinetic energy and the (r/r 0 ) N terms. We may therefore disregard the Coulomb, angular momentum and energy terms. Thus, we solve for the simplified differential equation

- 1 2
d 2 P (r) dr 2 + r r 0 N P (r) = 0, (32)
which for N = 2 is identical to the one-dimensional harmonic oscillator differential equation. Substituting

b = 2 r N 0 and c = N + 2 (33)
into Eq. ( 32) we get

d 2 P (r) dr 2 -b 2 r c-2 P (r) = 0. ( 34
)
With the transformation P (r) = √ rF (r) we can rewrite Eq. ( 34) as

r 2 d 2 dr 2 + r d dr -b 2 r c + 1 4 F (r) = 0. ( 35
)
Changing variables ξ = 2br c/2 /c and introducing α = 1/c, we obtain the modified Bessel's equation,

ξ 2 d 2 F (ξ) dξ 2 + ξ dF (ξ) dξ -ξ 2 + α 2 F (ξ) = 0. ( 36
)
Solutions are given in terms of the modified Bessel functions of the first and second kinds -that is, I α (ξ) and K α (ξ). Both are real when α is real and arg ξ = 0, which is the case here. Relevant asymptotic forms in our case are (see Eqs. 10.30.4, 10.25.3 in Ref. [START_REF] Dlmf | NIST Digital Library of Mathematical Functions[END_REF])

I α (z) ∼ e z √ 2πz , |arg z| < π 2 , (37)
K α (z) ∼ π 2z e -z , |arg z| < 3π 2 . ( 38
)
We may therefore exclude the modified Bessel function of the first kind since it is exponentially increasing.

After back substitution, we may write the asymptotic form for the physically acceptable solution of Eq. (32) as

P (r) ∼ r r 0 -N/4 exp - 2 √ 2r (N + 2) r r 0 N/2 . ( 39
)
For N = 0, we retrieve the expected Z-dependent exponential decay of the (free) hydrogenic wave function, whereas for N = 0, the Z-dependence vanishes and the asymptotics are fully driven by the confinement. For N = 2 (harmonic confinement), we obtain a Gaussian decay in agreement with the harmonic-oscillator wave function. For N > 2, we get stronger than Gaussian decay [48], suggesting that, for strong confinements, Gaussian basis sets are perhaps not so well suited for the problem and large basis-set expansions will be required. This will be analyzed in more detail below.

Figure 3 shows the 2s wave function for r 0 = 10a 0 and N = 10 in comparison with the respective functions of the free system and the hard-wall confinement. We see that spatial confinement pushes the wave function towards the nucleus, forcing it to approach zero more rapidly for r → +∞.

D. Convergence towards hard-wall results

Hard-wall confinement of the hydrogen atom has been reviewed by Aquino [29] and Laughlin [START_REF] Laughlin | [END_REF]. More recently, Katriel and Montgomery analyzed the convergence to the hard-wall limit for the confined hydrogen 16)). r0 = 10a0 (dashed line) and N = 10 is chosen.

For the soft-wall potential the asymptotic form of the wave function valid at large distances r defined in Eq. ( 39) is also shown.

atom at a given confinement radius [49]. Analytical solutions to the soft confinement problem is not known even for the hydrogen atom and we shall rather seek numerical solutions. However, placing the hard wall at the nodes of the free hydrogenic eigenstates gives eigenvalues and eigenstates for the confined system that are identical to analytical solutions of the free hydrogen atom [49]. It is not known how the position of the nodes of an eigenstate in a soft confinement potential depends on the stiffness parameter N . To analyze this in more detail, we discuss nodal features for the hard-wall potential first. Let P n : [0, ∞) → R be the solution of the radial Schrödinger equation of the free hydrogen atom in Eq. ( 16) with V c = 0. If n > + 1, then P n has n --

1 > 0 nodes at 0 < w n ,1 < w n ,2 < • • • < w n ,n--1 < ∞ where P n (w n ,i ) = 0. For each i ∈ {1, 2, . . . n --1}, the restriction P n | [0,w n ,i ] : [0, w n ,i ] → R of P n
from the non-negative axis [0, ∞) to the finite interval [0, w n ,i ] is the radial wave function of the hard-wall confined hydrogen atom with i -1 radial nodes and confinement radius r 0 = w n ,i . These confined systems have the same energy as the free hydrogen atom, E n = -1/(2n 2 ). Hence, for each state Ψ n m of the free hydrogen atom, there exist n--1 hard-wall confined systems with the same energy and the same wave function inside the walls. More generally, this result follows from a piecewise decomposition of the wave function with nodal boundaries as detailed by Courant and Hilbert [74]; see also Refs. [75,76].

We remind the reader that the integer n that appears in the discussion above is the principal quantum number of the free hydrogen atom but not of the confined atom. For all (confined or unconfined) atoms, the principal quantum number is equal to the total number of radial and angular nodes plus one. For a confined hydrogen atom of radius r 0 = w nl,i , the principal quantum number is thus given by n = i + and we get

E HW n , (w n ,i ) = E HW i+ , (w n ,i ) = E n . (40)
Note that, for the confined hydrogen atom, states of the same principal quantum number n but different angularmomentum quantum numbers are not degenerate, as discussed later in this section.

Letting n → ∞ in Eq. ( 40), using the continuity of E HW i+ , and the limit lim n→∞ E n = 0, we obtain

E HW i+ , (w ∞, ,i ) = 0. ( 41
)
For example, the critical radius at which the lowest energy of a hard-wall confined system of angular momentum becomes non-negative is therefore given by

r HW c, = w ∞, ,1 . (42)
The radial nodes at the n → ∞ limit can be obtained from the Sommerfeld-Welker relation [2,[START_REF] Ferreyra | [END_REF]77],

w ∞, ,i = 1 8 z 2 2 +1,i , (43)
where z 2 +1,i is the ith root of the Bessel function J 2 +1 (x) with x ≥ 0. With i = 1, from (42) we obtain for = 0, 1, 2, 3, respectively, the following increasing critical radii r HW c, in atomic units: 1.835246330, 5.088308227, 9.617366042, 15.36345002. Hence, r HW c,0 is the smallest nodal position that can be obtained for an analytical closed-form solution of the wave function. They do, however, arise from zeros of special confluent hypergeometric functions [START_REF] Scherbinin | Physics of Clusters[END_REF][START_REF] Ferreyra | [END_REF].

The question now arises as to how these nodes and the critical radius are influenced by soft confinement. To analyze numerically the convergence of the nodal positions, energies and critical radii in more general terms, we use Numerov's method [78,79] to obtain accurate solutions for Eq. ( 16) with V conf (r) = V N (r). This fourth-order linear multi-step method is used in standard applications to solve ordinary second-order differential equations. The method is implicit but can be made explicit if the differential equation is linear, as in our case. The efficiency of the Numerov method arises since we obtain a local error of O(h 6 ) with respect to step size h, with just one evaluation of the function per step. By comparison, the Runge-Kutta algorithm requires six function evaluations per step to achieve a local error of O(h 6 ) [START_REF] Süli | An Introduction to Numerical Analysis[END_REF].

The step size of the numerical calculations has been fixed to δr = 0.0001a 0 , whereas the upper integration limit for the range [0, r int ] was adjusted to the stiffness parameter chosen. For increasing stiffness parameter N , the integration limit r int was set closer to the confinement radius to avoid numerical instabilities. However, r int > r 0 was always kept as large as possible to ensure that the error due to the integration limit is below the numerical accuracy of Numerov's method. The hard-wall potential was obtained by setting the value of the last grid point to where E n (N, r SW c ) = 0 towards the hard-wall limit (bottom). The log-log inset plots showcase the power-law asymptotics.

zero. Our results are in perfect agreement with previous results by Aquino et al. [38].

Numerical calculations were also performed using a modified version of the finite-element-based atomic MCHF program LUCAS [START_REF] Sundholm | [END_REF]82]. The radial grid for hydrogen consisted of 200 non-equidistant elements, whose length increases exponentially with the radial distances, except for the first element. Within each element, the wave function is expanded using five equidistant Lagrange interpolation polynomials.

For a fixed radius r 0 and increasing stiffness parameter N , the nodal positions in the wave functions approach as expected the corresponding hard-wall nodes. Consider the difference between soft-and hard-wall nodes ∆w n ,i (r 0 , N ) = w SW n ,i (r 0 , N ) -w HW n ,i (r 0 ), (44) where w HW n ,i (r 0 ) = w SW n ,i (r 0 , N → ∞). In Figure 4, we illustrate the convergence of the soft-wall nodal positions towards the hard-wall nodes. It is clear that the convergence is slow -very large N values are required to approach the hard-wall limit. In general, the N -dependence of w SW n ,i is non-monotonic, becoming monotonic only with tighter confinement (smaller r 0 ). However, we note that, in the asymptotic region (N > 1000), the inequality w SW n ,i (r 0 , N ) > w HW n ,i (r 0 ) always holds, meaning that w SW n ,i converges towards the hard-wall limit from above. As is clearly seen from the log-log plots in Figure 4, the nodal position difference ∆w n ,i follows a simple powerlaw asymptotic behavior

∆w n ,i (r 0 , N ) = AN -B , (45)
where the exponent B typically has values in the range 0.8 -0.9. The exponent is virtually independent of quantum numbers n and (all asymptotic lines in Figure 4 are parallel) and shows a weak linear dependence on r 0

B(r 0 ) ≈ 0.87 -0.008r 0 . (46)
The deviations from the hard-wall results are much larger for smaller r 0 ; note the different scales in the plots of Figure 4.

The convergence of the energies to the hard-wall limit are shown in Figure 5. Just as for the nodal convergence, we observe generally non-monotonic behavior. Asymptotically, the soft-wall energy E SW n approaches E HW n from below and follows an inverse power-law

∆E n (r 0 , N ) = A -1 N B , (47)
where A and B are the same as in Eqs. ( 45) and ( 46) (for the same value of r 0 ). The magnitude of ∆E n increases as 1s < 2p < 2s < 3d < 3p < 3s, reflecting the increasing radial extents of the corresponding wave functions. For larger values of r 0 , the non-monotonic Ndependence stems from two competing effects. First, the confining potential V c influences the inner part of the wave function, raising its energy; with increasing N , this inner effect fades quickly and V c acts only in the vicinity of the wall. Second, the wave function is able to penetrate the soft wall, reducing its energy (relative to the hard-wall case); the higher the N , the less penetrable the wall becomes.

The soft-wall critical radius r SW c converges to r HW c from below, as illustrated in Figure 4. The curves corresponding to different states are shifted along the x axis (stiffness) and well separated by the quantum number n.

Figure 6 shows the convergence of the expectation values r -1 SW , r SW and r 2 SW towards their corresponding hard-wall limits. Whereas r -1 SW converges to the hard-wall limit from below, r SW and r 2 SW converge from above, as follows from the observation that the electron density is shifted towards the nucleus with increasing N .

For all shown properties, the asymptotic convergence towards the hard-wall limit is driven by a power function as in Eqs. ( 45) and (47). In all cases, the value of the exponent B stays fairly constant, whereas the prefactor A varies greatly between different properties and different states. The log-log inset plots showcase the power-law asymptotics.

Sen et al. have pointed out that the ith n-electron energy E i (Ω) obtained by solving the Schrödinger equation on the region Ω subject to the constraint that the wave function vanishes on the boundary δΩ (the Dirichlet boundary conditions) satisfies the relation [83].

Ω 1 Ω 2 =⇒ E i (Ω 1 ) > E i (Ω 2 ) . ( 48
)
This result applies to hard-wall confinement, where the wave function vanishes at the wall. For all quantum numbers, E HW n (r 0 ) therefore increases monotonically with decreasing r 0 . For soft confinement, the same result does not follow since we do not impose Dirichlet conditions on the wave function. However, for the soft-wall confinement, E n (r 0 ) > E n (r 0 + ) with > 0 holds to first order for each eigenstate of the hydrogen atom. Monotonicity should therefore hold in general if higher-order terms are small -that is, for large r 0 . This monotonic behavior is clearly seen in Figure 7, which also illustrates the slow convergence to the free hydrogen result for the harmonic confinement. Goldman and Joslin [34] discussed the dependence of the hydrogen energy levels on the angular-momentum quantum number for a fixed principal quantum number n in terms of the kinetic energy dominating the potential energy with decreasing r 0 . Their results are in agreement with Eq. ( 23), which holds for large N and therefore also in the hard-wall limit. Figure 7 illustrates the splitting of the 2s and 2p levels for the hard-wall potential and three different soft-wall potentials. We note the increasing separation of the two levels with decreasing r 0 . The energy ordering for higher principal quantum numbers is very similar -see, for example, Goldman and Joslin [34]. We note that an analysis of the nodal regions of the wave functions by Wilson implies E n, ≥ E n, -1 [76]. Hence, when the degeneracy of levels of the same principal quantum number n is lifted, then the state of higher angular momentum is lower in energy.

E. Using Gaussian Basis Sets

For hard-wall confinements, the multi-electron wave function must vanish at the wall and special basis functions must be designed to satisfy the Dirichlet boundary conditions [29,85]. This makes it difficult to study confinement effects using standard quantum-chemistry program packages. By not imposing such conditions, softwall confinement potentials are better suited for calculations with Gaussian basis functions. In fact, and most conveniently, the local pseudopotential operator [86] can be used by setting the Gaussian exponent to zero, which is much simpler than the implementation of the potential ansatz of Zicovich-Wilson [87]. On the other hand, since the wave function is pushed towards the nucleus and decays rapidly near r 0 , standard Gaussian basis sets may not be well suited for soft-wall confined systems with large N and small r 0 [48]. To make this point clear, we study the hydrogen atom using an even-tempered Gaussian basis. For chosen values of r 0 and N , two parameters α > 0 and β > 1 set up a geometric progression for the Gaussian exponents [88]:

γ k = αβ k , k = 1, . . . , M (49)
separately for each angular momentum . Even-tempered basis sets are asymptotically complete when α and β are allowed to depend on M [89][90][91].

We employ a non-linear conjugate gradient method [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF] to optimize α and β [START_REF] Research | Mathematica, Version 8.0[END_REF] using an energy threshold of < 10 -10 E h . However, with increasing N and decreasing r 0 , the threshold had to be lowered for numerical stability as the energy of the confined system becomes large. We ran calculations with r 0 = 0.5, . . . , 10a 0 and N ∈ {2, 5, 10, 20}. Different pairs of (α i , β i ) were explored through random selection in predefined intervals to ensure that the global minimum is obtained for a given confinement radius resulting in a final set of (α min (N, r 0 ), β min (N, r 0 )) values. The results using 20 Gaussian functions are displayed in Figure 8.

With increasing stiffness N , the optimized β min shifts to higher values, so that the even-tempered basis spans a broader range of Gaussian exponents γ. For small values of r 0 , the optimized α min increases with the increasing stiffness N effectively shifting the entire basis towards harder exponents. For higher r 0 , the trend is reversed and the basis shifts towards more diffuse exponents with increasing stiffness.

As seen in the bottom plot of Figure 8, for N = 2, the energy error introduced by the Gaussian basis is virtually independent of r 0 and is only governed by numerical noise, so that the Gausssian basis calculations perfectly reproduce the numerical results. This is not surprising considering the fact that the asymptotic behavior of the wave function is Gaussian for harmonic potentials. With increasing stiffness N , the basis set error increases and it grows considerably as r 0 approaches 0.

Our basis-set studies imply that small values of β and many basis functions with large exponents are required to describe properties over the whole range of confinement radii and stiffness parameters. For the dipole polarizability of hydrogen shown in Figure 9 calculated solving the Coupled-Perturbed Hartree-Fock equations [START_REF] Langhoff | [END_REF], we therefore used a much larger even-tempered 35s23p basis set (with center of basis functions at γ 18,s = 104.847 and γ 12,p = 2.643986 and fixed scaling factors β s = β p = 1.6) for stiffness parameters N = 2, 6 and 10. We note that variations of the γ values had very little effect on the results. The β p values could not be further reduced because of numerical problems for small radii r 0 .

As seen from Figure 9, the dependence of the polarizability on r 0 is similar for soft and hard confinements. For large radii, the curves approach the exact non-relativistic value of 4.5 a.u. of the free hydrogen atom; for small radii, all curves converge to zero polarizability. However, the confinement radius r 0 for which the dipole polarizability starts to change significantly depends critically on the stiffness parameter N . The harmonic potential (N = 2), in particular, exhibits a slow convergence towards the free hydrogenic value and shows a more gradual tran- sition between the r 0 → 0 and r 0 → ∞ limits. For larger values of r 0 , the stiffer confining potentials only act locally in the region where the wavefunction decays naturally with little to no effect on polarizability. On the other hand, the softer potentials (especially the harmonic potential) have significant effect in the region below the r 0 , thus lowering the polarizability even at large radii r 0 . Hence, the trend of polarizabilty versus stiffness reverses between the r 0 → 0 and r 0 → ∞ limits. More details on the hydrogenic results can be found in the supplementary Hard-wall results from Refs. [36,84]. All values are in atomic units.

material.

IV. SOFT-WALL CONFINED MANY-ELECTRON SYSTEMS

A. Computational details

We consider three many-electron atoms in our nonrelativistic confinement studies: helium ( 1 S), carbon (the three lowest states 3 P, 1 D and 1 S), and potassium (three lowest states 2 S, 2 P and 2 D). Methane is included as a molecule because of its high T d symmetry.

For the hard-wall confined helium, we performed numerical Hartree-Fock (HF) and multi-configuration HF (MCHF) calculations using a modified version of the finite-element atomic MCHF program LUCAS [START_REF] Sundholm | [END_REF]82] with the same grid as for the hydrogen atom. In the numerical complete active space self-consistent field (CASSCF) calculations the two electrons were distributed in 6s,5p,4d,3f,2g,1h active shells (here and in the following the numbers specify the number of shells, not quantum numbers). When the active space is increased further, the energies approach the values obtained by Montgomery Jr. et al. using explicitly correlated variational perturbation theory (VPT) [43].

For carbon and potassium, similar grids were used as for helium but with a radial grid consisting of 100 elements (401 grid points). For the carbon atom, the 1s orbital was inactive in our numerical CASSCF calculations with four electrons in 2s,2p,1d shells [denoted (1s//2s2p1d)]. In the numerical CASSCF calculations on potassium, the 12 electrons in the 1s-3s and 2p atomic orbitals were inactive and the remaining seven active electrons were distributed in 1s,2p,2d active shells [denoted (3s1p//1s2p2s)].

For the soft-wall confinement of helium, stiffness parameters N = 2, 6 and 10 were chosen. For the HF and CISD (FCI) calculations, special Gaussian type orbital (GTO) basis sets with many tight functions are required. For helium, even the uncontracted aug-cc-pV6Z correlation consistent basis set [95] is insufficient, the correlation energy plotted against r 0 showing unphysical wiggles. We therefore used a large uncontracted even-tempered basis set spanned by 33s23p17d13f11g9h functions with stiffness parameter N = 10. As for hydrogen, the energy is lowered with decreasing scaling factor β between the exponents. Because of linear dependencies, the smallest factor that could be used was 1.5. For each angular momentum, we chose α 0 small enough to be applicable for the stiffness parameters chosen (see supplementary material). Dipole polarizabilities were calculated numerically using the Finite-Field Perturbation Theory [96].

For the GTO calculations on potassium, we used an even-tempered basis with 41s27p19d functions. Stateaveraged CASSCF calculations with one electron in 1s,1p,1d active orbitals and all degenerate components of the lowest 2 S, 2 P, and 2 D states were included in the state averaging. This is equivalent to HF calculations for each state.

For CH 4 , the GTO basis-set requirements are somewhat less critical than for the atoms. We used the uncontracted aug-cc-pV6Z s and p sets [95] for C and H, augmented by the uncontracted aug-cc-pwCV5Z sets [97] for the higher angular-momentum functions. For the electron-correlation treatment, we used the coupledcluster singles-doubles-perturbative-triples (CCSD(T)) method. Unless otherwise noted, all electrons were correlated.

All GTO calculations were carried out using the MOL-PRO program package [START_REF] Werner | MOL-PRO, version 2018.1, A package of ab initio programs[END_REF].

B. The confined helium atom

The free helium isoelectronic series has been a test bench for electron-correlation treatments since the pioneering work of Hylleraas [69,70]. Upon application of the coordinate scaling r → r/Z mentioned in Section III B to the Hamiltonian of each member of the series and dividing by Z 2 , the Hamiltonian becomes

Ĥ(1, 2) = ĥZ=1 (1) + ĥZ=1 (2) + 1 Z 1 r 12 ( 50
)
where ĥZ=1 is the Hamiltonian of the hydrogen atom.

The factor λ = 1/Z appears as a natural perturbation parameter such that the energy of a two-electron atom of charge Z may be expressed as

E(Z) = Z 2 i=0 E (i) λ i = -Z 2 + 5 8 Z + E (2) + O(Z -1 )
where E (i) is the ith-order energy of the transformed Hamiltonian. The radius of convergence λ c = 1/Z c of the expansion corresponds to the critical nuclear charge Z c ∼ 0.91103 below which the atom is unstable [START_REF] Estienne | [END_REF]. The HF energy may be developed in an analogous manner [100]. An important observation is that the 1/Z expansions for the exact and HF energies are identical to first order, since the zeroth-order wave functions are the same [101]. Furthermore, since the second-order contribution to the energy in Eq. ( 50) is independent of Z, the correlation energy is independent of Z to lowest order and given by the limiting value

lim Z→∞ E corr (Z) = E (2) -E (2) HF . (51)
For the free atoms, the limiting value is E corr (Z → ∞) = -0.046663254E h [102]. These results may be extended to helium-like atoms under soft confinement provided that the confinement radius r 0 is scaled according to Eq. ( 31). This is demonstrated in Figure 10, where the CISD correlation energy of the helium isoelectronic series under soft confinement is plotted as a function of 1/Z, using a Z-dependent confinement radius r 0 (Z) = 2a 0 /Z 1+ 2 N . The curves are very nearly linear, as confirmed by linear regression, with slopes E

(3) corr = {0.0072706, 0.0070114, 0.0069044}E h for N = {2, 6, 10}. The corresponding intercepts are {-0.045967, -0.046841, -0.047207}E h and equal to the limiting correlation energies E

(2) corr . These energies may alternatively be calculated from an MP2-like expression [100,103] starting from the one-electron confined solutions, as we have confirmed by calculation.

We now examine the dependence of the (correlation) energy on the confinement parameters in more detail. In the top part of the Figure 11, we have plotted the CI ground-state electronic energy of helium against r 0 for different values of the stiffness parameter N . These energies are compared with the values obtained at the MCHF level using the hard-wall confinement potential. With increasing stiffness parameter, the CI energies slowly approach the energy of the hard-wall confined helium atom.

The CI energies are compared with energies calculated by Montgomery Jr. et al. using variational perturbation theory (VPT) [43] in the Supplementary Material.

For the correlation energy E corr = E CISD -E HF , studied for the hard-wall confinement by Gimarc in 1967 [104] and by Ludena and Gregory in 1979 [105], we observe an interesting behavior as shown in the lower panel of Figure 11. With increasing confinement radius r 0 , the correlation energy E corr (r 0 ) approaches the energy of the free helium atom -that is, E corr (r 0 → ∞) = -0.042044384E h for the exact value [106,107]. (The computed MCHF and CISD/GTO correlation energies are -0.041796E h , and -0.041919E h , respectively: the CISD value differs from the exact value because of basisset incompleteness, whereas the MCHF value differs because of the restricted orbital space.) Convergence to the free-helium correlation energy is much faster for larger stiffness parameters N since such potentials do not in-fluence the Coulomb potential significantly at distances r < r 0 .

Recalling that the high-density limit of correlation energy is finite, E corr (Z → ∞) = -0.046663254E h [102], we can expect to find a similar behavior when r 0 → 0. Indeed, Wilson et al. [108] showed computationally that the correlation energy in the hard-wall highdensity limit approaches that of ballium (two electrons confined in a sphere with potential V (r 12 ) = r -1 12 ) -that is, E ballium corr (r 0 → 0) = -0.055176E h [55], which can be compared with the correlation energy of -0.0545931E h obtained at the MCHF level for the hard-wall confined helium as r 0 → 0 using a 6s,5p,4d,3f,2g,1h active space. The good agreement between the correlation energy of ballium and helium when r 0 → 0 suggests that the correlation energy is independent of the nuclear charge in the limit of r 0 → 0. The discrepancy of -0.0005829E h between the correlation energy of ballium in the highdensity limit and the MCHF correlation energy in the limit r 0 → 0 is larger than the difference of -0.000371E h between the correlation energies calculated at the MCHF and VPT levels for r 0 = 1.0a 0 [43], suggesting that the convergence of the correlation energy is slightly slower for strongly confined atoms than for free atoms.

Analogously to the hard-wall case, the high-density limit of harmonically confined (N = 2) helium should approach that of 2-electron harmonium (Hookium) atom, i.e. two electrons confined by the harmonic potential, to which the Coulomb potential serves only as a perturbation. The limiting correlation energy value in this case is E Hookium corr (r 0 → 0) = -0.049703E h [109,110]. Harmonium systems with 2-6 electrons for both r 0 → 0 and r 0 → ∞ limits were previously extensively studied by Cioslowski and coworkers [110][111][112][113][114][115][116].

The soft-wall calculations show that the magnitude of the correlation energy in the r 0 → 0 limit increases with the increasing value of the stiffness parameter N , that is,

|E corr (N, r 0 → 0)| < |E corr (N + 1, r 0 → 0)|. ( 52
)
The same order of E corr (N, r 0 ) is apparent in the short range of the confinement radius r 0 , see Figure 11. This is consistent with the reported progression from harmonium to ballium [55,102]. The long-range behavior of the correlation energy including the maximum at about 2a 0 has been explained by Wilson et al. for hard-wall potentials to originate from different r 0 dependence of radial and angular correlation contributions [108]. The hard-wall numerical MCHF calculations show that the kinetic energy, the one-electron potential energy, and the two-electron repulsion energy for small r 0 are given by T (r 0 ) = 9.8832r -2 0 , V 1 (r 0 ) = -9.7481r -1 0 , and V 2 (r 0 ) = 1.7859r -1 0 , respectively. These scalings can be inferred from Eq.( 7) and translate into the well-known V -2/3 and V -1/3 scaling of kinetic and potential energy, respectively, with respect to volume (see for instance [117, sec. 15-4]).

Figure 12 shows the CI kinetic and potential energy contributions for helium as a function of the confinement radius r 0 for N = 10.

The r 2 expectation value and dipole polarizability at the CI level of theory are shown in Figure 13, where they are compared with the hard-wall values. They behave in a manner similar to that of the hydrogen atom. For a recent hard-wall treatment of the polarizability of raregas atoms, see Ref. [START_REF] Ludlow | [END_REF].

Experimental pressure-volume data for solid fcc/hcp helium are available from the work of Loubeyre et al. [119] and Mao et al. [120]. While the volume is given by V = (4π/3)r 3 0 for the hard-wall potential, it is not well-defined for a soft-wall potential since the wave function then extends to infinity. However, noting that the pressure for hard-wall confinement may be rewritten as

P c = - ∂E ∂V = - 1 4πr 2 0 ∂E ∂r 0 = E + T -V N 4πr 3 0 ( 53
)
where the last step follows from the virial theorem Eq. ( 11) [29], we shall take this expression as our definition of pressure in the case of soft-wall confinement.

If we consider a single atom under pressure, the P (V ) curves for the hard-or the soft-wall confinements will be markedly different from the real bulk system containing infinitely many atoms or molecules. In the latter case, the wave function of a single atom is allowed to spread out in the bulk material. Still, it is interesting to compare soft-wall results with the two limiting cases, the bulk and the hard-wall confined system. For this, we chose helium as an example. Hydrogen is more complicated because of strong covalent bonding at low pressures [121], but it has been considered by Ley-Koo and Rubinstein [9] and by LeSar and Herschbach [122]. The calculated P (V ) curves are shown in Figure 14, compared with the experimental data for solid fcc/hcp helium by Loubeyre et al. [119] and Mao et al. [120].

The P (V ) curves of the confined single atoms in Figure 14 are qualitatively similar to the corresponding solidstate curves, but are shifted to markedly higher pressures. In the long-range the soft van der Waals interaction dominates and solid He can be compressed with moderate pressures. At shorter He-He distances where overlap and repulsive forces dominate, further compression requires much higher pressure. Concerning the confined atom, the harmonic potential perturbs the Coulomb potential already at small distances, leading to high pressures in the long range. Increasing stiffness towards the hardwall confinement leads to smaller pressures in the long range but larger pressures in the short range. To simulate helium compression with confinement potentials will require higher stiffness parameters and a large number of helium atoms inside the confinement sphere.

C. The confined carbon atom

The energies for the three lowest-lying states of the confined carbon atom are shown in Figure 15. The ground-state configuration of the free carbon is Energy / E 1s 2 2s 2 2p 2 . Numerical CASSCF calculations on the three lowest states of C in a hard-wall potential show that confinement leads to changes in the occupation of the 2s and 2p shells; see Figure 15 inset. Since the 2s orbital with one radial node is larger than the 2p shell, for very small confinement radii the 1s 2 2p 4 configuration is significantly lower in energy than the 1s 2 2s 2 2p 2 configuration. It should be noted that the order of the 3 P, 1 D, and 1 S states is the same for all r 0 , because the three valence states have the same electron configuration for all confinement radii. Even though there are no state crossings, the configurations contributing to the wave functions change strongly. State crossings occur for potassium, because the electron configurations of the low-lying valence states differ, as discussed below.

h T -1/2 V 1 V N V 2 -1/2(V 1 +V 2 ) -1/2(V 1 +V 2 -NV N )

D. The confined potassium atom

The ground state of the free potassium atom is 2 S, the 2 P and 2 D states being the first and the second excited states, respectively. As shown in the Figure 16, at a soft-wall (N = 10) confinement radius of about 7a 0 the 2 P and 2 D states cross, the 2 D state becoming the first excited stated. The excitation energy to the 2 D state increases slightly upon compression until a radius of about 6a 0 , after which it decreases rapidly. At about 4.6a 0 , the 2 D state becomes the ground state, while the 2 S state rises steeply in energy, because the more compact nodeless 3d shell is less affected by compression than the more diffuse 4s orbital. Patil and Varshni found a similar crossing between 2 S and 2 P levels of Li within a hard-wall confinement using an effective Hamiltonian [71]. At small r 0 values, the kinetic energy dominates and states with rich radial nodal structure are energetically penalized. Further CASSCF calculations with larger active spaces, which included the 3p and 3d orbitals and more states, showed that at a confinement radius of about 2.6a 0 the 2 P state crosses further 2 P, 2 D, and 2 F states, which mainly arise from 3p 6 3d 1 4s 0 , 3p 5 3d 1 4s 1 , and 3p 5 3d 0 4s 2 configurations. At even smaller confinement radii, these states also cross the 2 S state. In this region, many states become nearly degenerate, and we were unable to unravel this complicated situation in detail. These calculations show that the electronic structure can significantly change when many-body atoms are exposed to high pressures. A similar situation is observed in strong magnetic fields, i.e., at field strengths of about one atomic unit B 0 = 235 kT.

E. The confined methane molecule

For CH 4 , the C-H optimized bond distance decreases monotonically with decreasing confinement radius r 0 , as shown in Figure 17. We only optimized the distances at the HF level, since the correlation effect is small relative to the effect of the confinement potential. The corresponding electronic energies (with all electrons correlated) as a function of the confinement radius r 0 are shown in Figure 18 (top).

The total correlation E corr and the core-core plus corevalence correlation E CC+CV corr energies for CH 4 (computed as the difference between the CCSD(T) energies with all electrons calculated and those with C(1s) uncorrelated) are shown in Figure 18 (bottom). As for helium, E corr goes through a maximum with decreasing r 0 , after which the absolute value steeply increases, most likely towards the high-density limit of the n-electron ballium system. As the molecule is compressed, also E CC+CV corr increases in magnitude.

In Figure 19, we show plots for the one-electron density of the free CH 4 molecule compared with that of the confined molecule with r 0 = 2.4a 0 and N = 10. We clearly see how the density becomes more spherical and compact under compression.

Finally, we considered the CH + 4 cation, using the same structures as for neutral CH 4 . A comparison between the total energies of CH + 4 and CH 4 shows that the two curves cross at a small confinement radius of about r 0 = 2.2a 0 (for N = 6, 10). To make this more apparent, we show the difference of the total energies of both molecules in Figure 20. The electron cannot escape from the confinement and remains bound despite the fact that CH + 4 is lower in energy. Nevertheless, we may regard this as the onset of a plasma-like behavior under such extreme conditions. This behavior is (perhaps) related to the fact that atomic or molecular solids can become conducting at extreme pressures. Similar effects have been found in confined atoms and were interpreted as ionization in confinement [123].

It should be noted, however, that the CH + 4 cation is Jahn-Teller distorted and has C 2v symmetry. This distortion lowers the energy of the cation, moving the crossing distance of the CH + 4 and CH 4 energies to larger r 0 values. Without the confinement potential, the H1-C-H2 and H3-C-H4 bond angles for CH + 4 are computed (CCSD(T)/V5Z with frozen core, HF values in parentheses) to be 125.65 • (124.9 • ) and 55.0 • (54.8 • ), and the respective CH bond distances are 1.081 Å (1.072 Å) and 1.185 Å (1.176 Å), respectively. With the confinement potential, the structure remains qualitatively the same, but for r 0 = 2.2a 0 , N = 10 the HF bond angles change to 117.9 • and 60.3 • and the distances are reduced to 0.768 Å and 0.794 Å, respectively (computed with the large basis set). It is interesting to note that H3-C-H4 then forms an almost equilateral triangle. Using the HF optimized structures for both species with r 0 = 2.2a 0 and N = 10, the CH + 4 ROHF-UCCSD(T) energy is 2.12 eV lower than the CH 4 one, while it is 0.98 eV higher if both are computed using the neutral CH 4 structure. Thus, at this confinement radius, the geometry relaxation lowers the CH + 4 energy by more than 3 eV.

V. CONCLUSIONS

We have undertaken detailed studies on the properties of the H, He, C and K atoms as well as the CH 4 and CH + 4 molecules confined within the spherical soft-wall potential

V N (r) = r r 0 N ( 54
)
of confinement radius r 0 and stiffness parameter N ; the hard-wall potential of radius r 0 is recovered as N tends to infinity. For hydrogen, we have carefully analyzed the use of Gaussian basis sets. While a harmonic potential is well suited for Gaussian functions, the expansion in Gaussian functions becomes increasingly more difficult with increasing stiffness parameter N . Having relatively wellbehaved basis sets for stiffness parameter N ≤ 10 avoiding near linear dependencies, we were able to present accurate shifts in the total energy and in the correlation energy relative to the free atom or molecule. We established a virial theorem for the soft confinements.

From our analysis, it is evident that the soft-wall potential is well-behaved for N and can easily be used within standard quantum-chemical program packages for many-electron atoms and molecules. Unlike the hardwall potential in Eq. (2), where an additional boundary condition is needed to force the wave function to become zero at the wall, the soft-wall potential results in standard one-electron integrals -already contained, for example, within matrix elements over the local pseudopotential operator [86]. The soft-wall confinement potential can easily be made finite by choosing the ansatz,

V N (r) = r r 0 N e -γr 2 , (55)
The soft-wall potential introduced in this work can be used, for example, in the simulation of materials under pressure. On the other hand, its extension to the relativistic domain is a challenge [46,124,[START_REF] Noon | A fully relativistic treatment of confined hydrogen-like atoms[END_REF]. These topics are the subject to our future studies.

See Supplemental Material at [URL will be inserted by publisher] for tables of properties for the confined systems studied.

FIG. 1 .

 1 FIG.1. The potential -1/r+VN (r) plotted with r0 = 1a0 and different stiffness parameters N . For a fixed r0, all potentials meet at (r0, 1 -1/r0) with slope 1/r 2 0 + N/r0.

FIG. 2 .

 2 FIG. 2. Expectation value r Nn of the free hydrogen atom plotted against the stiffness parameter N for n = 1, 2, 3. Red circles mark the -independent r -1 n values.

FIG. 3 .

 3 FIG.3. Radial part of the normalized 2s eigenfunctions of the free, soft-wall and hard-wall confined hydrogen atom (see Eq. (16)). r0 = 10a0 (dashed line) and N = 10 is chosen.For the soft-wall potential the asymptotic form of the wave function valid at large distances r defined in Eq. (39) is also shown.

FIG. 4 .

 4 FIG. 4. Convergence of the soft-wall node positions w SW n ,i towards the hard-wall limit for r0 = 10a0 (top) and r0 = 2a0 (center); convergence of the ground state critical radius r SW c

5 FIG. 5 .

 55 FIG.5. Convergence of soft-wall energy E SW towards the hard-wall limit for r0 = 10a0 (top) and r0 = 2a0 (bottom). The log-log inset plots showcase the power-law asymptotics.

FIG. 6 .

 6 FIG.6. Convergence of the soft-wall expectation values r -1 SW (top), r SW (center), r 2 SW (bottom) towards the hard-wall limit for r0 = 6a0. The log-log inset plots showcase the power-law asymptotics.

FIG. 7 .

 7 FIG.7. Lifting of the degeneracy of the hydrogen n = 2 levels by the confinement potential. The hard-wall data are from Refs.[34,37,60,84].

20 FIG. 8 .

 208 FIG. 8.Optimized even-tempered basis parameters αmin(r0, N ) (top) and βmin(r0, N ) (center) for a set of 20 Gaussian functions and the corresponding energy errors (bottom). Note that 20 basis functions are insufficient to describe high N / low r0 confinement.

FIG. 9 .

 9 FIG. 9. Dipole polarizability of the confined hydrogen atom with respect to the confinement radius r0 and stiffness N . Hard-wall results from Refs.[36,84]. All values are in atomic units.

FIG. 10 .

 10 FIG.10. CISD correlation energy of the confined helium isoelectronic series plotted as a function of 1/Z at r0(Z) = (2.0/Z 1+ 2 N )a0 for different values of the stiffness parameter N . In order to avoid basis set artifacts, the exponents of the 27s22p14d12f9g7h basis set for helium were scaled as αi(Z) = (Z/2) 2 αi(2).

FIG. 11 .

 11 FIG.11. Total ground-state CI energies (top) and the corresponding correlation energies (bottom) of the hard-and softwall confined helium atom with different stiffness parameters N . The hard-wall results are from numerical MCHF calculations with a 6s,5p,4d,3f,2g,1h active space.

FIG. 12 .

 12 FIG.12. Energy contributions for N = 10 (top); kinetic energy (second); N VN (third); and virial ratios -VC / 2T (bottom) for soft-wall confined helium. MCHF hard-wall data are included for comparison. hy.v. denotes the hypervirial ratio (-VC + N VN )/2 T derived from Eq. (9).

10 FIG. 13 .

 1013 FIG. 13. r 2 expectation value (top) and dipole polarizability (bottom) of the soft-wall confined helium atom for different stiffness parameters.

FIG. 14 .

 14 FIG.14. Pressure plotted against volume for the confined helium atom at the CI level for stiffness parameters N = 2, 6, 10 and at the MCHF hard-wall level, compared with measurements for the fcc and hcp phases by Loubeyre et al.[119] and Mao et al.[120].

a 0 FIG. 15 .

 015 FIG.15. Energies for the three lowest-lying states of carbon for a stiffness parameter N = 10 using a Gaussian basis. The inset shows the natural occupation numbers of the 2s and 2p shells for the 3 P state of carbon as a function of the hardwall confinement radius. The corresponding energies of the hard-wall confined carbon calculated at the CASSCF level are given in the Supplemental Material.

10 FIG. 17 .

 1017 FIG. 17. Hartree-Fock optimized C-H bond distance in CH4as a function of the confinement radius r0 for N = 2, 6, 10.

FIG. 18 .

 18 FIG. 18. Total CCSD(T) energy (top), total electron correlation energy Ecorr (bottom) and core-core plus core-valence correlation energy E CC+CV corr (bottom inset) for the soft-wall confined CH4 molecule plotted against the confinement radius r0 for N = 2, 6, 10. The symmetry was kept in T d .

FIG. 20 .

 20 FIG. 20. Difference in CCSD(T) energies of soft-wall confined CH + 4 and CH4 with changing confinement radius r0 for N = 2, 6, 10. Solid lines represent vertical ionization energies (CH4 geometry), dashed line shows adiabatic ionization energy (relaxed CH + 4 geometry).

ACKNOWLEDGMENTS

This work is in memory of Michael Wormit and Boris Pavlov, who both passed away at an early stage of this project in 2015. The authors thank Jonas Wiebke, Kyle Beloy, Dirk Andrae and Mathieu Lewin for fruitful discussions. PS is indebted to the Royal Society of New Zealand for financial support in terms of a Marsden Fund (14-MAU-034). DS acknowledges financial support from Magnus Ehrnrooth Foundation, Swedish Cultural Foundation in Finland, and the Academy of Finland (275845 and 314821). TH acknowledges support by the Norwegian Research Council through the CoE Hylleraas Centre for Quantum Molecular Sciences Grant No. 262695. LFP is grateful for the support from the Slovak Research and Development Agency (APVV-15-0105) and the Scientific Grant Agency of the Slovak Republic (1/0777/19). The authors are grateful to the Centre for Advanced Study at the Norwegian Academy of Science and Letters, Oslo, Norway, where much of this work was carried out under the project "Molecules in Extreme Environments" during the academic year 2017-2018.

Appendix A: Proof of r N n,l < r N n,l-1

Here, we present a proof by induction of the inequality r N n,l < r N n,l-1 , using the Kramers-Pasternack recurrence relation in Eq. (19). We use the induction step N -1 → N and consider the two cases 0 < N < 2l + 1 and N ≥ 2l + 1.

1. Case 1: 0 < N < 2l + 1

We start by dividing Eq. ( 19) with r N -1 n,l :

with

and proceed with induction for the ratios R

n,l-1 holds. For N < 2l + 1 the coefficient 1 4 N (2l+1+N )(2l+1-N ) is always positive and from Eq. (A1) we get

Further, we note that

and we have

Using Eq. (A1) we get

The inequality still holds if we drop the right-hand term, yielding

that concludes this part of the induction. Finally, we must make the connection with r N n,l < r N n,l-1 . From Eq. (A7), we follow

Again, we prove by induction, only this time directly for r N n,l < r N n,l-1 . For the initial step we know already from Case 1 that the inequality holds for N = 2l-1 and N = 2l. We assume that r N -2 n,l < r N -2 n,l-1 and r N -1 n,l < r N -1 n,l-1 . As for N ≥ 2l + 1, the rightmost term in Eq. ( 20) is non-negative and replacing r N -1 n,l and r N -2 n,l yields

Using Eq. (A4), we get

and using Eq. ( 19) gives

As in the last paragraph, the last term is always negative so we have

This concludes our proof.

Here, we prove the inequality

From Z-scaling it follows that

Hence, it suffices to prove

For the following, we set Z = 1. We apply Hölder's inequality

with p -1 + q -1 = 1 for the domain A ⊂ R d . Next, we introduce

and we get (the integrands are all positive functions)

dr

To prove that this is indeed the case, we start from the relation [62]

in terms of the Hahn polynomials

The generalized hypergeometric series is given by [72, eq.16. or in terms of gamma functions [72, eq.5.2.5]

In our case, we have r = s + 1 for the Hahn polynomials (B8). Due to appearance of -k < 0 in (B8), the power series will reduce to a polynomial of maximal order k.

Using the recurrence relation [72, 5.5.1]

we can rewrite the Hahn polynomials as

We simplify further using the relation

and get

We further have

We note that

as well as

Note also that r N +1 n has one more term in the expansion than r N n for Z = 1. From this, we conclude that

which completes our proof since r 0 n = 1. We finally mention that the more general inequality (B23)