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