pairinteraction
A Rydberg Interaction Calculator
SystemAtom.cpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2024 Pairinteraction Developers
2// SPDX-License-Identifier: LGPL-3.0-or-later
3
5
12
13#include <Eigen/Dense>
14#include <algorithm>
15#include <limits>
16#include <memory>
17#include <set>
18#include <spdlog/spdlog.h>
19
20namespace pairinteraction {
21template <typename Scalar>
22SystemAtom<Scalar>::SystemAtom(std::shared_ptr<const basis_t> basis)
23 : System<SystemAtom<Scalar>>(std::move(basis)) {}
24
25template <typename Scalar>
26SystemAtom<Scalar> &SystemAtom<Scalar>::set_electric_field(const std::array<real_t, 3> &field) {
27 this->hamiltonian_requires_construction = true;
28
29 if (!traits::NumTraits<Scalar>::is_complex_v && field[1] != 0) {
30 throw std::invalid_argument(
31 "The field must not have a y-component if the scalar type is real.");
32 }
33
34 electric_field = field;
35
36 return *this;
37}
38
39template <typename Scalar>
40SystemAtom<Scalar> &SystemAtom<Scalar>::set_magnetic_field(const std::array<real_t, 3> &field) {
41 this->hamiltonian_requires_construction = true;
42
43 if (!traits::NumTraits<Scalar>::is_complex_v && field[1] != 0) {
44 throw std::invalid_argument(
45 "The field must not have a y-component if the scalar type is real.");
46 }
47
48 magnetic_field = field;
49
50 return *this;
51}
52
53template <typename Scalar>
55 this->hamiltonian_requires_construction = true;
56 diamagnetism_enabled = enable;
57 return *this;
58}
59
60template <typename Scalar>
62SystemAtom<Scalar>::set_ion_distance_vector(const std::array<real_t, 3> &vector) {
63 this->hamiltonian_requires_construction = true;
64
65 if (!traits::NumTraits<Scalar>::is_complex_v && vector[1] != 0) {
66 throw std::invalid_argument(
67 "The distance vector must not have a y-component if the scalar type is real.");
68 }
69
70 ion_distance_vector = vector;
71
72 return *this;
73}
74
75template <typename Scalar>
77 this->hamiltonian_requires_construction = true;
78 ion_charge = charge;
79 return *this;
80}
81
82template <typename Scalar>
84 this->hamiltonian_requires_construction = true;
85 if (value < 2 || value > 3) {
86 throw std::invalid_argument("The order of the Rydberg-ion interaction must be 2 or 3");
87 }
88 ion_interaction_order = value;
89 return *this;
90}
91
92template <typename Scalar>
94 auto basis = this->hamiltonian->get_basis();
95
96 // Construct the unperturbed Hamiltonian
97 this->hamiltonian = std::make_unique<OperatorAtom<Scalar>>(basis, OperatorType::ENERGY);
98 this->hamiltonian_is_diagonal = true;
99 bool sort_by_quantum_number_f = true;
100 bool sort_by_quantum_number_m = true;
101 bool sort_by_parity = true;
102
103 // Estimate the numerical precision so that we can decide which terms to keep
104 Eigen::VectorX<real_t> diag = this->hamiltonian->get_matrix().diagonal().real();
105 real_t scale = (diag - diag.mean() * Eigen::VectorX<real_t>::Ones(diag.size())).norm();
106 real_t numerical_precision = 100 * scale * std::numeric_limits<real_t>::epsilon();
107
108 real_t typical_magnetic_dipole = 1e2; // ~n^1
109 real_t typical_electric_dipole = 1e4; // ~n^2
110 real_t typical_electric_quadrupole = 1e8; // ~n^4
111
112 // Add the interaction with the field of an ion (see
113 // https://en.wikipedia.org/wiki/Multipole_expansion#Spherical_form for details)
114
115 Eigen::Map<const Eigen::Vector3<real_t>> vector_map(ion_distance_vector.data(),
116 ion_distance_vector.size());
117 real_t distance = vector_map.norm();
118 SPDLOG_DEBUG("Distance to the ion: {}", distance);
119
120 // Dipole order
121 if (std::isfinite(distance) && ion_interaction_order >= 2) {
122 // Calculate sqrt(4pi/3) * r * Y_1,q with q = -1, 0, 1 for the first three elements. Take
123 // the conjugate of the result and scale it by 1/distance^2.
124 Eigen::Vector3<Scalar> vector_dipole_order =
125 spherical::get_transformator<Scalar>(1) * vector_map / distance;
126 vector_dipole_order = vector_dipole_order.conjugate() / std::pow(distance, 2);
127
128 for (int q = -1; q <= 1; ++q) {
129 if (std::abs(vector_dipole_order[q + 1]) * typical_electric_dipole >
130 numerical_precision) {
131 *this->hamiltonian -= ion_charge * vector_dipole_order[q + 1] *
133 this->hamiltonian_is_diagonal = false;
134 sort_by_quantum_number_f = false;
135 sort_by_parity = false;
136 sort_by_quantum_number_m &= (q == 0);
137 }
138 }
139 }
140
141 // Quadrupole order (the last entry of vector_quadrupole_order would correspond to
142 // ELECTRIC_QUADRUPOLE_ZERO and does not contribute to a *traceless* quadrupole)
143 if (std::isfinite(distance) && ion_interaction_order >= 3) {
144 // Calculate sqrt(4pi/5) * r^2 * Y_2,q / 3 with q = -2, -1, 0, 1, 2 for the first five
145 // elements and (x^2 + y^2 + z^2) / 6 as the last element. Take the conjugate of the
146 // result and scale it by 1/distance^3.
147 Eigen::Vector<Scalar, 6> vector_quadrupole_order = spherical::get_transformator<Scalar>(2) *
148 Eigen::KroneckerProduct(vector_map / distance, vector_map / distance);
149 vector_quadrupole_order = 3 * vector_quadrupole_order.conjugate() / std::pow(distance, 3);
150
151 for (int q = -2; q <= 2; ++q) {
152 if (std::abs(vector_quadrupole_order[q + 2]) * typical_electric_quadrupole >
153 numerical_precision) {
154 *this->hamiltonian -= ion_charge * vector_quadrupole_order[q + 2] *
155 OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, q);
156 this->hamiltonian_is_diagonal = false;
157 sort_by_quantum_number_f = false;
158 sort_by_quantum_number_m &= (q == 0);
159 }
160 }
161 }
162
163 // Add external fields (see https://arxiv.org/abs/1612.08053 for details)
164
165 // Transform the electric field to spherical coordinates and take the conjugate
166 Eigen::Vector3<Scalar> electric_field_spherical = spherical::get_transformator<Scalar>(1) *
167 Eigen::Map<const Eigen::Vector3<real_t>>(electric_field.data(), electric_field.size());
168 electric_field_spherical = electric_field_spherical.conjugate();
169
170 // Transform the magnetic field to spherical coordinates and take the conjugate
171 Eigen::Vector3<Scalar> magnetic_field_spherical = spherical::get_transformator<Scalar>(1) *
172 Eigen::Map<const Eigen::Vector3<real_t>>(magnetic_field.data(), magnetic_field.size());
173 magnetic_field_spherical = magnetic_field_spherical.conjugate();
174
175 // Stark effect: - \vec{d} \vec{E} = - d_{1,0} E_{0} + d_{1,1} E_{-} + d_{1,-1} E_{+}
176 // = - d_{1,0} E_{0} - d_{1,1} E_{+}^* - d_{1,-1} E_{-}^*
177 // with the electric dipole operator: d_{1,q} = - e r sqrt{4 pi / 3} Y_{1,q}(\theta, \phi)
178 // where electric_field_spherical=[E_{-}^*, E_{0}, E_{+}^*]
179 for (int q = -1; q <= 1; ++q) {
180 if (std::abs(electric_field_spherical[q + 1]) * typical_electric_dipole >
181 numerical_precision) {
182 *this->hamiltonian -= electric_field_spherical[q + 1] *
183 OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_DIPOLE, q);
184 this->hamiltonian_is_diagonal = false;
185 sort_by_quantum_number_f = false;
186 sort_by_parity = false;
187 sort_by_quantum_number_m &= (q == 0);
188 }
189 }
190
191 // Zeeman effect: - \vec{\mu} \vec{B} = - \mu_{1,0} B_{0} + \mu_{1,1} B_{-} + \mu_{1,-1} B_{+}
192 // = - \mu_{1,0} B_{0} - \mu_{1,1} B_{+}^* - \mu_{1,-1} B_{-}^*
193 // with the magnetic dipole operator: \vec{\mu} = - \mu_B / \hbar (g_l \vec{l} + g_s \vec{s})
194 // where magnetic_field_spherical=[B_{-}^*, B_{0}, B_{+}^*]
195 for (int q = -1; q <= 1; ++q) {
196 if (std::abs(magnetic_field_spherical[q + 1]) * typical_magnetic_dipole >
197 numerical_precision) {
198 *this->hamiltonian -= magnetic_field_spherical[q + 1] *
199 OperatorAtom<Scalar>(basis, OperatorType::MAGNETIC_DIPOLE, q);
200 this->hamiltonian_is_diagonal = false;
201 sort_by_quantum_number_f = false;
202 sort_by_quantum_number_m &= (q == 0);
203 }
204 }
205
206 // Diamagnetism: 1 / (8 m_e) abs(\vec{d} \times \vec{B})^2
207 // = (1/12) \left[ B_0^2 (q0 - d_{2,0}) + B_+ B_- (-2 q0 - d_{2,0})
208 // + \sqrt{3} B_0 B_- d_{2,1} + \sqrt{3} B_0 B_+ d_{2,-1}
209 // - \sqrt{3/2} B_-^2 d_{2,2} - \sqrt{3/2} B_+^2 d_{2,-2} \right]
210 // with the operator: q0 = e^2 r^2 sqrt{4 pi} Y_{0,0} = e^2 r^2
211 // and the electric quadrupole operator: d_{2,q} = e^2 r^2 sqrt{4 pi / 5} Y_{2,q}(\theta, \phi)
212 // where magnetic_field_spherical=[B_{-}^*, B_{0}, B_{+}^*]
213 if (diamagnetism_enabled) {
214 if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
215 numerical_precision) {
216 *this->hamiltonian += static_cast<real_t>(1 / 12.) *
217 static_cast<Scalar>(std::pow(magnetic_field_spherical[1], 2)) *
218 OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, 0);
219 *this->hamiltonian -= static_cast<real_t>(1 / 12.) *
220 static_cast<Scalar>(std::pow(magnetic_field_spherical[1], 2)) *
221 OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, 0);
222 this->hamiltonian_is_diagonal = false;
223 sort_by_quantum_number_f = false;
224 }
225 if (std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
226 numerical_precision &&
227 std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
228 numerical_precision) {
229 *this->hamiltonian -= static_cast<real_t>(2 / 12.) * magnetic_field_spherical[0] *
230 magnetic_field_spherical[2] *
231 OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, 0);
232 *this->hamiltonian -= static_cast<real_t>(1 / 12.) * magnetic_field_spherical[0] *
233 magnetic_field_spherical[2] *
234 OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, 0);
235 this->hamiltonian_is_diagonal = false;
236 sort_by_quantum_number_f = false;
237 }
238 if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
239 numerical_precision &&
240 std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
241 numerical_precision) {
242 *this->hamiltonian -= static_cast<real_t>(std::sqrt(3.0) / 12.) *
243 magnetic_field_spherical[1] * magnetic_field_spherical[2] *
244 OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, 1);
245 this->hamiltonian_is_diagonal = false;
246 sort_by_quantum_number_f = false;
247 sort_by_quantum_number_m = false;
248 }
249 if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
250 numerical_precision &&
251 std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
252 numerical_precision) {
253 *this->hamiltonian -= static_cast<real_t>(std::sqrt(3.0) / 12.) *
254 magnetic_field_spherical[1] * magnetic_field_spherical[0] *
255 OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, -1);
256 this->hamiltonian_is_diagonal = false;
257 sort_by_quantum_number_f = false;
258 sort_by_quantum_number_m = false;
259 }
260 if (std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
261 numerical_precision) {
262 *this->hamiltonian -= static_cast<real_t>(std::sqrt(1.5) / 12.) *
263 static_cast<Scalar>(std::pow(magnetic_field_spherical[2], 2)) *
264 OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, 2);
265 this->hamiltonian_is_diagonal = false;
266 sort_by_quantum_number_f = false;
267 sort_by_quantum_number_m = false;
268 }
269 if (std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
270 numerical_precision) {
271 *this->hamiltonian -= static_cast<real_t>(std::sqrt(1.5) / 12.) *
272 static_cast<Scalar>(std::pow(magnetic_field_spherical[0], 2)) *
273 OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, -2);
274 this->hamiltonian_is_diagonal = false;
275 sort_by_quantum_number_f = false;
276 sort_by_quantum_number_m = false;
277 }
278 }
279
280 // Store which labels can be used to block-diagonalize the Hamiltonian
281 this->blockdiagonalizing_labels.clear();
282 if (sort_by_quantum_number_f) {
283 this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_F);
284 }
285 if (sort_by_quantum_number_m) {
286 this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_M);
287 }
288 if (sort_by_parity) {
289 this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_PARITY);
290 }
291}
292
293// Explicit instantiations
294template class SystemAtom<double>;
295template class SystemAtom<std::complex<double>>;
296} // namespace pairinteraction
Type & set_ion_charge(real_t charge)
Definition: SystemAtom.cpp:76
typename traits::CrtpTraits< Type >::real_t real_t
Definition: SystemAtom.hpp:42
Type & set_diamagnetism_enabled(bool enable)
Definition: SystemAtom.cpp:54
SystemAtom(std::shared_ptr< const basis_t > basis)
Definition: SystemAtom.cpp:22
Type & set_electric_field(const std::array< real_t, 3 > &field)
Definition: SystemAtom.cpp:26
Type & set_magnetic_field(const std::array< real_t, 3 > &field)
Definition: SystemAtom.cpp:40
Type & set_ion_distance_vector(const std::array< real_t, 3 > &vector)
Definition: SystemAtom.cpp:62
Type & set_ion_interaction_order(int value)
Definition: SystemAtom.cpp:83
std::shared_ptr< const basis_t > get_basis() const
Definition: System.cpp:76
Matrix< Type, Size, 1 > Vector
Matrix< Type, 3, 1 > Vector3
Matrix< Type, Dynamic, 1 > VectorX
Helper struct to extract types from a numerical type.
Definition: traits.hpp:35