April 17, 2019 at 02:35

Muon catalysis in terms of quantum chemistry. Part II: electronic vs. muon chemical bond

    Mnogabukaff that quantum chemistry thinks about the principle of muon catalysis: how exactly muon lowers the temperature of the desired plasma. In two parts (the first part can be read here ).

    The essence of the second part is simple: the muon is heavier than the electron, so it provides a stronger chemical bond and a closer approach of the nuclei, thereby lowering the required plasma temperature for igniting the thermonuclear reaction.

    But those who want to look at the formulas, graphs, and see the conceptual essence of quantum chemistry as applied to the simplest (quasi) molecules, welcome under cat.


    Introduction


    In the first part (see here ) we examined the difference between a hydrogen atom$ \ mathrm {H} ^ \ cdot = \ mathrm {p ^ + e ^ -} $ from its heavy muon counterpart $ \ mathrm {p} ^ + \ mu ^ - $: in the second case, the muon will be tied more strongly, and it will sit at a closer distance from the proton. At the same time, we examined some important things that we will need here (forms of orbitals and the atomic system of units).

    In the second part (i.e., here) we will try to understand why, how and how much the plasma temperature required to ignite the thermonuclear reaction will decrease. The reactions that interest us are:

    $ \ mathrm {{} ^ nH + {} ^ m H \ rightarrow} \ \ text {new kernels + energy} $


    where n, m = 1,2,3 correspond to proton, deuterium and tritium, respectively. Naturally, these nuclei have a positive charge, so if you try to bring them closer, they will begin to repel in accordance with the Coulomb law (see the previous part ), and this is the very barrier that prevents the onset of fusion reactions. By the way, in the case of nuclear decay reactions, this repulsion has the opposite role, because after separation from the common nucleus, the fragments, repelling from each other, acquire additional kinetic energy, and it is this energy that is heated in nuclear power plants.

    To overcome this Coulomb barrier, an increase in plasma temperature ( T ) is required , which, as everyone remembers from the MKT school courseis related to the average particle velocity in plasma ( v ) by the formula

    $ mv ^ 2 = 3k_ \ mathrm {B} T $

    where m is the mass of particles, and$ k_ \ mathrm {B} $- Boltzmann constant .

    But, let’s imagine that we have combined two hydrogen nuclei into a certain particle, where they are already located close, and therefore the rest of the barrier for them is already very small. Then we would need to significantly accelerate these particles (read: we need lower temperatures) in order to combine them into something new. And just such a role should play an intermediate ion$ (\ mathrm {{} ^ nH} \ mu ^ - \ mathrm {{} ^ mH}) ^ + $, an analogue of the ion of a hydrogen molecule $ \ mathrm {H} _2 ^ + = (\ mathrm {He ^ -H}) ^ + $.
    Having examined the differences between these two particles, we will realize how effective the muon is in lowering the ignition temperature of thermonuclear fusion.

    MILK MO LKAO method


    So, we have our molecular system, consisting of 2 hydrogen nuclei with a charge + e (one electron charge modulo) and one particle (electron or muon) with a charge - e . Our system, until it collides with other particles, is isolated, and therefore its energy can be decomposed into its constituent parts:

    $ E = T (\ mathrm {H} _1) + T (\ mathrm {H} _2) + \ underbrace {T (\ mathrm {e} ^ - / \ mu ^ -) + V (\ mathrm {H_1 \ text {from} H_2}) + V (\ mathrm {\ mathrm {e} ^ - / \ mu ^ - \ text {к} \ mathrm {H} _1}) + V (\ mathrm {\ mathrm {e} ^ - / \ mu ^ - \ text {k} \ mathrm {H} _2})} _ {E_ \ mathrm {e}} $


    where the first two terms ($ T (\ mathrm {H} _1) $ and $ T (\ mathrm {H} _2) $) Is the kinetic energy of hydrogen nuclei, the third term ($ T (\ mathrm {e} ^ - / \ mu ^ -) $) Is the kinetic energy of a negative particle (electron or muon), the fourth term $ V (\ mathrm {H_1 \ text {from} H_2}) $Is the energy of the Coulomb repulsion of hydrogens from each other, and the remaining two are the Coulomb attraction of the electron / muon to each of the protons. In the general case, this is a 3-body problem, just a quantum one. Naturally, solving it in the forehead is very difficult. But, fortunately, the nuclei are at least 1800 times heavier than the electron, and 10 heavier than the muon, so they will move clearly slower than small negative particles. Due to this, you can first solve the problem in turn: first, find the energy of motions that are not related to the movement of nuclei, i.e.$ E_ \ mathrm {e} $and then full energy. It looks like this.

    1. The arrangement of hydrogen nuclei relative to each other is selected, and this determines the Coulomb interactions between them and with the electron / muon. Coulomb potential$ V (R) = k \ frac {q_1q_2} {R} $ depends only on particle charges $ q_i $and the distance between them, so for all hydrogen isotopes this value will be the same. Further, the problem of the motion of an electron / muon in the field of these nuclei is solved. This is the task of one body.
    2. These energies $ E_ \ mathrm {e} $are calculated for all possible arrangements of the nuclei relative to each other, and this will be the effective potential energy of motion of the nuclei. In our case, we need to calculate the energies at different distances relative to each other, so the potential for a pair of nuclei is always one-dimensional. Well, then we only need to solve the two-body problem of the motion of two hydrogen isotopes relative to each other.

    Obviously, the root of the problem with us is the calculation of the electron / muon energy in the field of nuclei $ E_ \ mathrm {e} $. In fact, this is the chemical bond: a certain potential that holds the nuclei together in certain places. And this very task of finding the energy of chemical bonding is the main one in quantum chemistry.

    Unfortunately, both the muon and the electron are quantum particles, therefore, in order to find this energy, we have to resort to the methods of quantum mechanics. In fact, our problem of the motion of an electron / muon in the field of two identical nuclei is solved explicitly (see here ), but this solution is very complicated and the result is not as clear as in the case of a hydrogen-like atom. Therefore, we will try to disassemble a different, approximate approach, which is applicable to any systems. This is the so-called molecular orbitals method as linear combinations of atomic orbitals, or MO LKAO.

    Let's take a closer look at the Schrödinger equation for the motion of an electron / muon in the field of hydrogen nuclei:

    $ \ hat {H} \ psi = \ underbrace {\ left (\ overbrace {- \ frac {1} {2m} (\ frac {\ partial ^ 2} {\ partial x ^ 2} + \ frac {\ partial ^ 2} {\ partial y ^ 2} + \ frac {\ partial ^ 2} {\ partial z ^ 2})} ^ {\ hat {T}} + \ overbrace {- \ frac {1} {R_1}} ^ {\ hat {V} _1} + \ overbrace {- \ frac {1} {R_2}} ^ {\ hat {V} _2} + \ overbrace {\ frac {1} {R}} ^ {\ hat {V } _ \ mathrm {HH}} \ right)} _ {\ hat {H}} \ psi = E \ psi $


    This equation was written in the atomic system of units (see PS in the previous part ), therefore the charge of the hydrogen nucleus and the electron / muon is +1, --1, respectively, the electron mass is m = 1, and for the muon m ≈207.

    And if you take a closer look, you can see that in the Hamiltonian you can select a piece connected purely with the movement of a negative particle around only one of the nuclei, which is just the Hamiltonian of the hydrogen atom, and this can be done in 2 ways:

    $ \ hat {H} = (\ overbrace {\ hat {T} + \ hat {V} _1} ^ {\ hat {H} _ {1}} + \ hat {V} _2 + \ hat {V} _ \ mathrm {HH}) = (\ overbrace {\ hat {T} + \ hat {V} _2} ^ {\ hat {H} _ {2}} + \ hat {V} _1 + \ hat {V} _ \ mathrm {HH}) $


    Outside the Hamiltonian of a hydrogen-like atom ($ \ hat {H} _i, \ i = $ 1.2) we always have 2 pieces: the energy of interaction of an electron / muon with another nucleus ($ \ hat {V} _j $) and nuclear repulsion energy ($ \ hat {V} _ \ mathrm {HH} $) The second of them does not affect the movement of electrons at all - it's just a shift of energy by a certain amount, but the interaction of an electron with another nucleus is an important thing.
    We can imagine that at any moment our particle rotates only around one of the nuclei, and the interaction with the second is just a correction. As a method of rotation around one of the nuclei, we can assume that the electron / muon is in the ground (1s) state, the wave function for which is well known to us from the previous part:

    $ | 1s \ rangle = \ frac {1} {\ sqrt {\ pi}} \ exp \ left (- \ frac {R} {R_1} \ right) $


    Where $ R_1 $Is the Bohr radius for a particle. In the case of an electron$ R_1 = 1 $ Boron (which is the Bohr radius for an electron, equal to about 0.5 angstroms), and in the case of a muon $ R_1 = \ frac {1} {m_ \ mu} \ approx \ frac {1} {207} $.
    In order to somehow approximate the electron / muon wave function in the field of 2 nuclei, we can try to take the following representation:

    $ \ psi \ approx c_1 | 1s_1 \ rangle + c_2 | 1s_2 \ rangle $


    and then the problem of solving a complex partial differential equation with us is reduced to searching for 2 unknown coefficients c 1 and c 2 . This is the very molecular orbital presented as the sum with the coefficients (a linear combination of the scientific) atomic 1s orbitals.

    Naturally, we need an equation for these parameters. And getting it is quite simple if you substitute this approximation in the Schrödinger equation$ \ hat {H} \ psi = E \ psi $:

    $ \ hat {H} (c_1 | 1s_1 \ rangle + c_2 | 1s_2 \ rangle) = c_1 \ hat {H} | 1s_1 \ rangle + c_2 \ hat {H} | 1s_2 \ rangle = E (c_1 | 1s_1 \ rangle + c_2 | 1s_2 \ rangle) = c_1 E | 1s_1 \ rangle + c_2 E | 1s_2 \ rangle $


    Actually, we want this ratio to be satisfied everywhere, so we can somehow calculate the average values ​​of this all. We multiply this equation on the left by$ <1s_1 | $ and $ <1s_2 | $and integrate over all coordinates. As a result, we obtain a system of 2 linear equations, where it is necessary to find the coefficients c 1 , c 2 and energy E :

    $ \ begin {pmatrix} \ langle 1s_1 | \ hat {H} | 1s_1 \ rangle & \ langle 1s_1 | \ hat {H} | 1s_2 \ rangle \\ \ langle 1s_2 | \ hat {H} | 1s_1 \ rangle & \ langle 1s_2 | \ hat {H} | 1s_2 \ rangle \\ \ end {pmatrix} \ begin {pmatrix} c_1 \\ c_2 \ end {pmatrix} = E \ begin {pmatrix} \ langle 1s_1 | 1s_1 \ rangle & \ langle 1s_1 | 1s_2 \ rangle \\ \ langle 1s_2 | 1s_1 \ rangle & \ langle 1s_2 | 1s_2 \ rangle \\ \ end {pmatrix} \ begin {pmatrix} c_1 \\ c_2 \ end {pmatrix} $


    Anyone who has studied linear algebra will recognize a generalized eigenvector-eigenvalue problem. Before solving it, we will analyze what the elements of the 2 matrices, the existing tuta, are equal to (and at the same time we introduce their short designation with one letter).

    • Let's start with the simplest: $ \ langle 1s_1 | 1s_1 \ rangle = \ langle 1s_2 | 1s_2 \ rangle = 1 $ - this is the normalization of wave functions, and as we recall, the total probability of finding an electron / muon is at least 1 somewhere.
    • $ \ langle 1s_1 | 1s_2 \ rangle = \ langle 1s_2 | 1s_1 \ rangle = S $- this is the so-called overlap integral, showing how 1s electron clouds overlap for each of the atoms.
    • $ \ langle 1s_1 | \ hat {H} | 1s_1 \ rangle = \ langle 1s_2 | \ hat {H} | 1s_2 \ rangle = \ alpha $. This integral consists of several parts:

      $ \ langle 1s_1 | \ hat {H} | 1s_1 \ rangle = \ underbrace {\ langle 1s_1 | \ hat {H} _1 | 1s_1 \ rangle} _ {- \ frac {m} {2}} + \ langle 1s_1 | \ hat {V} _2 | 1s_1 \ rangle + \ frac {1} {R} $


    • $ \ langle 1s_1 | \ hat {H} | 1s_2 \ rangle = \ langle 1s_2 | \ hat {H} | 1s_1 \ rangle = \ beta $. Here it is similar:

      $ \ langle 1s_2 | \ hat {H} | 1s_1 \ rangle = \ underbrace {\ langle 1s_2 | \ overbrace {\ hat {H} _1 | 1s_1 \ rangle} ^ {- \ frac {m} {2} | 1s_1 \ rangle}} {{\ frac {m} {2} S} + \ langle 1s_2 | \ hat {V} _2 | 1s_1 \ rangle + \ frac {S} {R} $


      those. the energy of a hydrogen-like atom and internuclear repulsion, scaled by the overlap integral (first and last terms), and, as it were, the energy of electron / muon hopping from one atom to another.

    Let us find the expressions for the energies of our hydrogen-like ion from the equation rewritten as

    $ \ begin {pmatrix} \ alpha & \ beta \\ \ beta & \ alpha \ end {pmatrix} \ begin {pmatrix} s_1 \\ c_2 \ end {pmatrix} = E \ begin {pmatrix} 1 & S \\ S & 1 \ end {pmatrix} \ begin {pmatrix} c_1 \\ c_2 \ end {pmatrix} $


    To find the energy you need to solve the equation:

    $ \ det \ begin {pmatrix} \ alpha -E & \ beta -ES \\ \ beta - ES & \ alpha -E \ end {pmatrix} = (\ alpha -E) ^ 2 - (\ beta - ES) ^ 2 = 0 $


    where "det" denotes the determinant (determinant of a matrix, in Russian).

    The solutions of this quadratic equation with respect to E are:

    $ E_ \ pm = \ frac {\ alpha \ pm \ beta} {1 \ pm S} = - \ frac {m} {2} + \ frac {1} {R} + \ frac {\ langle 1s_1 | \ hat {V} _2 | 1s_1 \ rangle \ pm \ langle 1s_2 | \ hat {V} _2 | 1s_1 \ rangle} {1 \ pm S} $


    The first piece is obviously the energy of the atom, the second is the internuclear repulsion, the same Coulomb barrier that prevents the ignition of the thermonuclear reaction, and the last complex structure must be dealt with.

    If we discard the internuclear repulsion, which is just a reference point for the electron / muon energy, we get that we have two states with energy

    $ \ epsilon_ \ pm = - \ frac {m} {2} + \ frac {\ langle 1s_1 | \ hat {V} _2 | 1s_1 \ rangle \ pm \ langle 1s_2 | \ hat {V} _2 | 1s_1 \ rangle} {1 \ pm S} $


    Since both wave functions $ | 1s_1 \ rangle $ and $ | 1s_2 \ rangle $ - positive, and $ \ hat {V} _i <0 $ (because the negative particle is always drawn to the positive), then $ \ epsilon_ + <- \ frac {m} {2} $ (energy of a single atom), and $ \ epsilon_-> - \ frac {m} {2} $, i.e. we get a standard picture of molecular orbitals:

    image

    Lower orbital with energy$ E _ + $ called binding, and the top (with energy $ E _- $) - anti-binding, or loosening. As a result, if an electron / muon sits on the lower molecular orbital, then it benefits from flying around 2 nuclei than around one, and with its movement it lowers the total energy of the system. And this is the very magical chemical bond that screens the internuclear repulsion, allowing the nuclei to be next to each other for quite some time.

    And here the integrals of the chemical bond must be calculated in order to understand how closely the hydrogen nuclei are allowed to be. In fact, all three of the sought-after integrals are calculated analytically, but it is terribly hemorrhoidal and complicated (anyone interested, see Chapter 9 in Flary’s Quantum Chemistry book ). Therefore, we will go a different way, simpler, and calculate these integrals numerically using the Monte Carlo method.

    Metropolis Method


    Мне видится очень логичным, в тексте про термоядерную энергию воздать почести её дедушке: военному атому, а если конкретнее, Манхэттенскому проекту. Именно из него вырос метод Монте-Карло, и в частности алгоритм Метрополиса, один из авторов которого, Эдвард Теллер, «отец водородной бомбы» (то бишь человек, запустивший термоядерный синтез на атолле Эниветок).

    Про принцип Арнольда
    Despite the fact that the algorithm bears the name of Metropolis, the last (and main type) author of the original article “Equation of State Calculations by Fast Computing Machines” in the Journal of Chemical Physics , various sources claim that the author of the algorithm is Enrico Fermi, rediscovered by Stan Ulam (both are not in the authors at all), or Marshal Rosenbluth, who is in third place in the list of authors, while Nicholas Metropolis is the author of the term "Monte Carlo" and the person who provided computer time . All in full accordance with the principle of Arnold . : D

    In general, we will analyze the essence of the method. It is intended for the tasks of statistical mechanics. The main distribution in it is the Boltzmann distribution: the probability of detecting a system in a certain state is$ \ exp (- \ beta E) $, $ \ beta ^ {- 1} = k_ \ mathrm {B} T $. And the observed value of some parameter A for the system in thermodynamic equilibrium is equal to the integral

    $ \ langle A \ rangle = \ frac {1} {Z} \ int A (q) \ exp (- \ beta E (q)) dq $


    where q is the coordinates that parameterize the state of the system (for example, the coordinates / momenta of the particles), and Z is the normalization factor called the partition function:

    $ Z = \ int \ exp (- \ beta E (q)) dq $


    If there are sooooo many particles in the system, then counting none of the integrals in the forehead is completely unrealistic. The naive Monte Carlo method, in which we simply select a bunch of random q coordinates , will also not give anything meaningful if there are really possible states of the system for which the probability$ \ exp (- \ beta E) $noticeably non-zero, very few. And it is precisely for such cases that we need a sample by significance, in which we allow the algorithm to sample only sufficiently probable places in the state space.

    The Metropolis algorithm looks as follows.

    1. When initiating the simulation, we select some starting approximation in the configuration space $ \ mathbf {q} ^ {(0)} $ and some vector of the maximum possible increment $ \ delta \ mathbf {q} $. At the starting point, we calculate the energy of the system$ E ^ {(0)} = E (\ mathbf {q} ^ {(0)}) $ (read - probability $ p = \ exp (- \ beta E ^ {(0)}) $)
    2. The new configuration at the nth step is as follows.
      1. Calculate the energy of the trial configuration $ E_ \ mathrm {trial} = E (\ mathbf {q} _ \ mathrm {trial}) $ (i.e. probability $ p_ \ mathrm {trial} = \ exp (- \ beta E_ \ mathrm {trial}) $)
      2. And then we compare the old probability $ p ^ {(n)} $ with trial $ p_ \ mathrm {trial} $

        • if the new configuration has a greater or the same probability ($ \ frac {p_ \ mathrm {trial}} {p ^ {(n)}} \ geq 1 $), or, equivalently, the energy of the new point is lower or the same as in the old ($ E_ \ mathrm {trial} \ leq E ^ {(n)} $), then the new point is accepted and the system goes into it ($ q ^ {(n + 1)} = q_ \ mathrm {trial} $),
        • if the trial configuration is higher in energy ($ E_ \ mathrm {trial}> E ^ {(n)} $), which is equivalent $ \ frac {p_ \ mathrm {trial}} {p ^ {(n)}} <1 $, then in this case we generate a random number $ P \ in [0; 1) $from a uniform distribution, and compare it with the ratio of the probabilities, which are the transition probabilities. If$ P <\ frac {p_ \ mathrm {trial}} {p ^ {(n)}} $, then we accept a new point, and if not ($P \geq \frac{p_\mathrm{trial}}{p^{(n)}} $), then we reject, and the system remains in the old configuration ($q^{(n+1)} = q^{(n)}$) ...


    3. Taking many steps according to the algorithm above, we sample a significant (i.e., really important) part of the possible space of system configurations. The integral of interest to us is calculated by the formula:
      $\langle A \rangle = \frac{1}{Z} \int A(\mathbf{q}) \exp(-\beta E(\mathbf{q})) d\mathbf{q} = \frac{1}{N}\sum_{n=0}^{N} A(\mathbf{q}^{(n)}) $

    This is how the Metropolis algorithm works.

    And now it would be necessary to adapt it to the calculation of the 3 integrals that interest us. Let's look at them in more detail.

    • $ S(R) = \langle 1s_2 |1s_1\rangle = \int \limits_{-\infty}^{+\infty} \int \limits_{-\infty}^{+\infty} \int \limits_{-\infty}^{+\infty} \overbrace{ \underbrace{ \frac{1}{\sqrt{\pi}} \exp(-m\underbrace{|\mathbf{r} - \mathbf{r}_2|}_{R_2}) }_{1s_2} }^{A(\mathbf{r})} \cdot \overbrace{ \underbrace{ \frac{1}{\sqrt{\pi}} \exp(-m\underbrace{|\mathbf{r} - \mathbf{r}_1|}_{R_1}) }_{1s_1} }^{p(\mathbf{r})} d x dy dz $where $\mathbf{r}=(x,y,z)^\mathbf{T}$ - coordinates of the electron / muon, $\mathbf{r}_i=(x_i,y_i,z_i)^\mathbf{T}$ Are the coordinates of the hydrogen nuclei, and $R_i = |\mathbf{r} - \mathbf{r}_i| = \sqrt{(x-x_i)^2 + (y-y_i)^2 + (z-z_i)^2} $ - distances between positive and negative particles,
    • $\langle 1s_1 |\hat{V}_2 |1s_1 \rangle = - \int \limits_{-\infty}^{+\infty} \int \limits_{-\infty}^{+\infty} \int \limits_{-\infty}^{+\infty} \overbrace{ \underbrace{ \frac{1}{\sqrt{\pi}} \exp(-mR_1) }_{1s_1} \frac{1}{R_2} }^{A(\mathbf{r})} \cdot \overbrace{ \underbrace{ \frac{1}{\sqrt{\pi}} \exp(-mR_1) }_{1s_1} }^{p(\mathbf{r})} dxdydz $
    • $\langle 1s_2 |\hat{V}_2 |1s_1 \rangle = - \int \limits_{-\infty}^{+\infty} \int \limits_{-\infty}^{+\infty} \int \limits_{-\infty}^{+\infty} \overbrace{ \underbrace{ \frac{1}{\sqrt{\pi}} \exp(-mR_2) }_{1s_2} \frac{1}{R_1} }^{A(\mathbf{r})} \cdot \overbrace{ \underbrace{ \frac{1}{\sqrt{\pi}} \exp(-mR_1) }_{1s_1} }^{p(\mathbf{r})} dxdydz $

    It can be seen that, if we calculate the 1s-function of one of the atoms for the probability p ,

    to do so, of course, is not very good,
    because the probability density is the modulus of the squared wave function $|\psi|^2$, not the wave function itself $\psi$.

    then everything else under the sign of the integral (the second wave function and in 2 out of 3 cases the potential of attraction of the electron / muon to the nucleus) will be a function whose average value is calculated. The only thing that will have to be done, in contrast to the usual calculation by the Metropolis method, is to straighten the normalization of the integrals. The fact is that the standard normalization will be on

    $Z=\int \limits_{-\infty}^{+\infty} \int \limits_{-\infty}^{+\infty} \int \limits_{-\infty}^{+\infty} \exp(-m R) dx dy dz = 4 \pi \int \limits_{0}^{+\infty} \exp(-m R) R^2 dR = \frac{8 \pi}{m^3}$


    And we need normalization to $\sqrt{\langle 1s_1 | 1s_1 \rangle}$where

    $\langle 1s_1 | 1s_1 \rangle = \int \limits_{-\infty}^{+\infty} \int \limits_{-\infty}^{+\infty} \int \limits_{-\infty}^{+\infty} \exp(-2 m R) dx dy dz = 4 \pi \int \limits_{0}^{+\infty} \exp(-2m R) R^2 dR = \frac{\pi}{m^3}$


    This means that each integral calculated according to the Metropolis will need to be multiplied by a factor

    $\frac{Z}{\sqrt{\langle 1s_1 | 1s_1 \rangle} } = 8 \sqrt{\frac{\pi}{m^3}}$



    This can already be arranged in the form of a certain script, for example, in Python (for example, the code is below).

    For example, like that.
    import numpy as np
from math import *
# r     = 0...+infty
# phi   = 0...2pi
# theta = 0...pi
# function to convert spherical coordinates into Cartesian
def sph2cart(r, phi, theta):
    xyz=np.zeros(3)
    xyz[0]=r*cos(phi)*sin(theta)
    xyz[1]=r*sin(phi)*sin(theta)
    xyz[2]=r*cos(theta)
    return xyz
# Distance between vectors r1 and r2
def dist(r1, r2):
    return sqrt(np.dot(r1-r2, r1-r2))
# re -- Cartesian coordinates of electron
# rn -- Cartesian coordinates of nucleus
# psi_1s returns a value of 1s wavefunction
def psi_1s(re, rn, scale=1.0):
    ren=dist(re,rn)
    return scale**(3/2)*exp(-scale*ren)/(sqrt(pi))
########################################
############ Settings ##################
########################################
NumPtsPerIntegral=100000   # Number of points per integral... duuuh
#mass=1.0                     # mass of the particle (electron)
mass=207.0                  # mass of the particle (muon)
AllRab=[]
#AllRab+=[1.0*(0.1)**n for n in range(0,10) ]
AllRab+=[ (1.4+0.25*n)/mass for n in range(0,10)]
print(AllRab)
########################################
########################################
########################################
dumpster=open("res.dat", "w") # output file to store results of the simulation
dx=2.0/mass   # maximal increment for the coordinate
rna=np.array([0.0, 0.0, 0.0])     # position of nucleus "a"
renorm=8.0*sqrt(pi/mass**3)  # factor to readjust result from incorrect norm of Metropolis weighting to a correct 1s wavefunction norm
# loop for the potential energy calculation at the chosen internuclear distances
for npt,Rab in enumerate(AllRab): 
    Norm=0.0  # <1s_a | 1s_a > for check
    Sab=0.0   # <1s_a | 1s_b >
    Vaa=0.0   # <1s_a | |r - R_b|**(-1) | 1s_a >
    Vab=0.0   # <1s_a | |r - R_b|**(-1) | 1s_b >
    re=np.array([1.0/mass, 0.0, 0.0]) # initial position of the electron
    rnb=np.array([Rab, 0.0, 0.0])     # position of nucleus "b"
    NumAcc=0.0                        # Number of accepted points
    for i in range(0,NumPtsPerIntegral):  # loop for Metropolis algorithm 
        newre=re+np.random.uniform(low=-dx, high=dx, size=3)  # trial position of electron
        pnew=psi_1s(newre, rna, scale=mass) ## trial probability
        pold=psi_1s(re, rna, scale=mass)    ## previous probability due to dumb and ineffective realization
        if pnew/pold >= np.random.random(): ## importance sampling step 
            re=newre
            NumAcc+=1. 
        Norm+=psi_1s(re, rna, scale=mass)   
        Sab+=psi_1s(re, rnb, scale=mass)    
        Vaa+=psi_1s(re, np.zeros(3), scale=mass)/dist(re, rnb)
        Vab+=psi_1s(re, rnb, scale=mass)/dist(re, rnb)
    Norm*=renorm/NumPtsPerIntegral
    Sab*=renorm/NumPtsPerIntegral
    Vaa*=renorm/NumPtsPerIntegral
    Vab*=renorm/NumPtsPerIntegral
    def s_test(x,scale=1.0):  # this is an analytical expression for overlap integral S in case of 1s hydrogen wavefunctions
        return exp(-x*scale)*(1.+x*scale+(1./3.)*(x*scale)**2)  
    #E=-0.5*mass**2 + 1./sqrt(np.dot(rnb,rnb)) - (Vaa + Vab)/(1.0 + Sab)       # full energy
    E= 1./sqrt(np.dot(rnb,rnb)) - (Vaa + Vab)/(1.0 + Sab)    # energy adjusted to energy of a single atom as the dissociational limit
    dumpster.write(" %10.3e %40.10f %15.10f %15.10f   %15.5f  %15.5f\n"  % (Rab, E, Sab, s_test(Rab,scale=mass), 100.0*NumAcc/NumPtsPerIntegral, Norm )) 
    dumpster.flush()
dumpster.close()



    Using such calculations, we can finally compare the potential energies in the hydrogen ion $\mathrm{H}_2^+$ and its muon counterpart.

    $\mathrm{H_2^+ = p^+e^- p^+}$ vs. $\mathrm{p^+ \mu^- p^+}$


    So, armed with a script, we can calculate the surface of the potential energy of the approach of the hydrogen nuclei bound by an electron and a muon. As an energy reference point, we take atoms infinitely diluted from each other (i.e.,$-m/2$, which is equal to the potential at the distance between the nuclei $R = + \infty$)

    In the case of an electron, the potential near the minimum looks like this: The



    minimum occurs at a distance of about 2 Boron (i.e., approximately the sum of 2 atomic radii), and the energy of dissociation of the molecule into fragments is approximately 0.06 Hartree, which corresponds to heating to about 20,000 degrees Kelvin (or Celsius) , it doesn’t matter). To convert energies, I recommend using online resources, such as this .

    A similar situation with the muon bound hydrogen ion:



    Since the Bohr radius for muon hydrogen is smaller (see the previous part), then the hydrogen nuclei in the minimum potential energy sit approximately 200 times closer. The breakdown energy of this molecule is already more than 10 Hartree, which corresponds to a temperature of more than three lyam degrees ($\approx (3.2 \cdot 10^6 )^\circ$)

    For ignition, the reactions usually require a temperature of the order of 10 8 K , which is about 320 Hartree. Let us see at what distances a similar energy is achieved in the case of the ordinary divodoron ion and in the case of its muon version:



    In the first case, this corresponds to a distance of about 0.0058 Boron (vertical line).

    A similar distance in muon hydrogen is achieved at an energy of about 190 Ha, i.e. about one and a half times less. And this is the simplest estimate of the temperature of muon catalysis.

    But in fact, everything will be even cooler. The fact is that if a stable particle is formed$\mathrm{({}^mH (\mu^-) {}^nH)^+}$, then these nuclei, while the muon is alive, will oscillate relative to each other. And here tunneling from the “two hydrogen atoms” state to the “heavier core” state can occur, and the probability of tunneling depends on the required tunneling length d approximately as$p^{-d}$, so that bringing the two nuclei closer together by the muon we will greatly increase the probability of the tunneling course of this reaction. Unfortunately, estimates of this effect no longer require quantum chemistry, but nuclear physics, so this part of the discussion is beyond the scope of this post. So on this we will stop.

    PS Why is it not so simple?


    In fact, to form these particles is not so simple in plasma conditions. The fact is that if we collide two particles, then their total energy obviously exceeds the energy of dissociation (or ionization, in the case of a nucleus + electron / muon), so when they collide, they do not form a stable particle (atom, ion, molecule), but fly by past each other. To stick to each other, they need to throw off somewhere surplus energy, and for this we need a third extra who will take on this energy. It can be a photon, or some kind of left particle flying nearby, but the main thing is that the conditions should contribute to this entrainment of excess energy.

    PPS


    If you have any comments / clarifications / questions, write in comments or in PM. I will correct everything, I will answer and explain everything.

    Only registered users can participate in the survey. Please come in.

    How complicated is the material presented or its presentation?

    • 9% Afftar kill yourself up the wall !!! 6
    • 31.8% White little fluffy fox !!! How difficult!!! 21
    • 24.2% Oh, I remember with a squeak my student years ... 16
    • 19.6% It seems you can live. thirteen
    • 30.3% The most, write more! 20
    • 15.1% Kindergarten! : D 10
    66 users voted. 18 users abstained.
    Tags:

    Also popular now: