Line data Source code
1 : # SPDX-FileCopyrightText: 2025 PairInteraction Developers
2 : # SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : """Tests for Rydberg ions, i.e., charged species whose interaction includes monopole terms."""
5 :
6 1 : from __future__ import annotations
7 :
8 1 : from typing import TYPE_CHECKING
9 :
10 1 : import numpy as np
11 1 : import pytest
12 1 : from pairinteraction.units import ureg
13 :
14 : if TYPE_CHECKING:
15 : from .utils import PairinteractionModule
16 :
17 :
18 1 : def test_monopole_matrix_elements(pi_module: PairinteractionModule) -> None:
19 : """The monopole operator is the total charge in units of the charge -e of the Rydberg electron."""
20 1 : ket_ion = pi_module.KetAtom("Sr88_ion", n=60, l=0, j=0.5, m=0.5)
21 1 : assert ket_ion.get_matrix_element(ket_ion, "electric_monopole", q=0, unit="e") == pytest.approx(-1)
22 :
23 : # The monopole operator is diagonal and the same for all states of the ion
24 1 : basis_ion = pi_module.BasisAtom("Sr88_ion", n=(59, 61), l=(0, 2))
25 1 : monopole = basis_ion.get_matrix_elements(basis_ion, "electric_monopole", q=0, unit="e").toarray()
26 1 : np.testing.assert_allclose(monopole, -np.eye(basis_ion.number_of_states))
27 :
28 :
29 1 : @pytest.mark.parametrize("order", [2, 3])
30 1 : @pytest.mark.parametrize("direction", ["z", "x", "y"])
31 1 : def test_atom_ion_pair_vs_point_charge(
32 : pi_module: PairinteractionModule, use_real: bool, order: int, direction: str
33 : ) -> None:
34 : """A Rydberg ion in a single state must act on a Rydberg atom like a classical point charge."""
35 1 : if direction == "y" and use_real:
36 1 : pytest.skip("a y-component requires complex numbers")
37 :
38 1 : distance = 3 # micrometer
39 1 : distance_vector = {"x": [distance, 0, 0], "y": [0, distance, 0], "z": [0, 0, distance]}[direction]
40 :
41 1 : basis_atom = pi_module.BasisAtom("Rb", n=(58, 62), l=(0, 3), m=(0.5, 0.5) if direction == "z" else None)
42 1 : ket_ion = pi_module.KetAtom("Sr88_ion", n=60, l=0, j=0.5, m=0.5)
43 1 : basis_ion = pi_module.BasisAtom("Sr88_ion", n=(0, 0), additional_kets=[ket_ion])
44 1 : assert basis_ion.number_of_states == 1
45 :
46 : # Pair system of the atom and the ion
47 1 : basis_pair = pi_module.BasisPair([pi_module.SystemAtom(basis_atom), pi_module.SystemAtom(basis_ion)])
48 1 : assert basis_pair.number_of_states == basis_atom.number_of_states
49 1 : system_pair = (
50 : pi_module.SystemPair(basis_pair)
51 : .set_interaction_order(order)
52 : .set_distance_vector(distance_vector, unit="micrometer")
53 : .diagonalize(diagonalizer="eigen")
54 : )
55 1 : energies_pair = system_pair.get_eigenenergies(unit="GHz") - ket_ion.get_energy(unit="GHz")
56 :
57 : # Atom in the field of a classical point charge
58 1 : system_reference = (
59 : pi_module.SystemAtom(basis_atom)
60 : .set_ion_charge(1, unit="e")
61 : .set_ion_interaction_order(order)
62 : .set_ion_distance_vector(distance_vector, unit="micrometer")
63 : .diagonalize(diagonalizer="eigen")
64 : )
65 1 : energies_reference = system_reference.get_eigenenergies(unit="GHz")
66 :
67 1 : np.testing.assert_allclose(energies_pair, energies_reference, atol=1e-6, rtol=0)
68 :
69 : # The interaction must have a significant effect
70 1 : energies_unperturbed = pi_module.SystemAtom(basis_atom).diagonalize(diagonalizer="eigen").get_eigenenergies("GHz")
71 1 : assert np.linalg.norm(energies_pair - energies_unperturbed) > 1e-3
72 :
73 :
74 1 : def test_ion_ion_coulomb_repulsion(pi_module: PairinteractionModule) -> None:
75 : """Two ions in a single state each interact via the repulsive Coulomb interaction Z1*Z2/R."""
76 1 : ket_ion = pi_module.KetAtom("Sr88_ion", n=60, l=0, j=0.5, m=0.5)
77 1 : basis_ion = pi_module.BasisAtom("Sr88_ion", n=(0, 0), additional_kets=[ket_ion])
78 1 : system_ion = pi_module.SystemAtom(basis_ion)
79 1 : basis_pair = pi_module.BasisPair([system_ion, system_ion])
80 1 : assert basis_pair.number_of_states == 1
81 :
82 1 : distances = np.array([0.5, 1, 2, 5]) # micrometer
83 1 : for order in [1, 2, 3]:
84 1 : energies = np.array(
85 : [
86 : pi_module.SystemPair(basis_pair)
87 : .set_interaction_order(order)
88 : .set_distance(d, unit="micrometer")
89 : .get_hamiltonian(unit="hartree")[0, 0]
90 : for d in distances
91 : ]
92 : )
93 1 : energies -= 2 * ket_ion.get_energy(unit="hartree")
94 1 : distances_au = ureg.Quantity(distances, "micrometer").to("bohr").magnitude
95 1 : np.testing.assert_allclose(energies, 1 / distances_au, rtol=1e-10)
96 1 : assert np.all(energies > 0)
97 :
98 :
99 1 : def test_ion_ion_multipole_orders(pi_module: PairinteractionModule) -> None:
100 : """The contributions of the individual multipole orders scale as 1/R^order."""
101 1 : ket = pi_module.KetAtom("Sr88_ion", n=60, l=0, j=0.5, m=0.5)
102 1 : basis = pi_module.BasisAtom("Sr88_ion", n=(ket.n - 1, ket.n + 1), l=(0, 2))
103 1 : system = pi_module.SystemAtom(basis)
104 :
105 1 : delta_energy = 300 # GHz
106 1 : pair_energy = 2 * ket.get_energy(unit="GHz")
107 1 : basis_pair = pi_module.BasisPair(
108 : [system, system], energy=(pair_energy - delta_energy, pair_energy + delta_energy), energy_unit="GHz", m=(1, 1)
109 : )
110 1 : assert basis_pair.number_of_states > 1
111 :
112 1 : distances = np.linspace(1, 5, 5)
113 1 : hamiltonian_0 = pi_module.SystemPair(basis_pair).get_hamiltonian(unit="GHz").toarray()
114 1 : hamiltonians = {0: [hamiltonian_0 for _ in distances]}
115 1 : for order in [1, 2, 3]:
116 1 : hamiltonians[order] = [
117 : pi_module.SystemPair(basis_pair)
118 : .set_interaction_order(order)
119 : .set_distance(d, unit="micrometer")
120 : .get_hamiltonian(unit="GHz")
121 : .toarray()
122 : for d in distances
123 : ]
124 1 : contributions = {
125 : order: np.linalg.norm(np.array(hamiltonians[order]) - np.array(hamiltonians[order - 1]), axis=(1, 2))
126 : for order in [1, 2, 3]
127 : }
128 :
129 1 : for order, norm in contributions.items():
130 1 : assert norm[0] > 0
131 1 : np.testing.assert_allclose(norm * distances**order, norm[0] * distances[0] ** order, rtol=1e-8)
132 :
133 : # The monopole-monopole interaction is a constant shift of all pair states
134 1 : shift = hamiltonians[1][0] - hamiltonians[0][0]
135 1 : np.testing.assert_allclose(shift, shift[0, 0] * np.eye(basis_pair.number_of_states), atol=1e-10)
136 :
137 :
138 1 : def test_neutral_pair_unaffected(pi_module: PairinteractionModule) -> None:
139 : """The monopole terms vanish for neutral atoms and must not alter the block structure of the Hamiltonian."""
140 1 : basis = pi_module.BasisAtom("Rb", n=(59, 61), l=(0, 2))
141 1 : system = pi_module.SystemAtom(basis)
142 1 : ket = pi_module.KetAtom("Rb", n=60, l=0, j=0.5, m=0.5)
143 1 : pair_energy = 2 * ket.get_energy(unit="GHz")
144 1 : basis_pair = pi_module.BasisPair([system, system], energy=(pair_energy - 3, pair_energy + 3), energy_unit="GHz")
145 :
146 : # For a distance vector along z, the Hamiltonian must be block-diagonal with respect to the quantum number m
147 1 : system_pair = pi_module.SystemPair(basis_pair).set_interaction_order(3).set_distance(3, unit="micrometer")
148 1 : basis_cpp = system_pair._cpp.get_basis()
149 1 : m = np.array([basis_cpp.get_quantum_number("m", i) for i in range(basis_cpp.get_number_of_states())])
150 1 : assert len(np.unique(m)) > 1
151 1 : hamiltonian = system_pair.get_hamiltonian(unit="GHz").toarray()
152 1 : assert np.all(hamiltonian[m[:, None] != m[None, :]] == 0)
153 :
154 : # Without charged species, the interaction orders 1 and 2 do not contribute
155 1 : system_pair_order_2 = pi_module.SystemPair(basis_pair).set_interaction_order(2).set_distance(3, unit="micrometer")
156 1 : system_pair_order_0 = pi_module.SystemPair(basis_pair)
157 1 : hamiltonian_order_2 = system_pair_order_2.get_hamiltonian(unit="GHz").toarray()
158 1 : hamiltonian_order_0 = system_pair_order_0.get_hamiltonian(unit="GHz").toarray()
159 1 : np.testing.assert_allclose(hamiltonian_order_2, hamiltonian_order_0)
|
