Line data Source code
1 : // SPDX-FileCopyrightText: 2024 Pairinteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/system/SystemAtom.hpp"
5 :
6 : #include "pairinteraction/enums/OperatorType.hpp"
7 : #include "pairinteraction/enums/TransformationType.hpp"
8 : #include "pairinteraction/operator/OperatorAtom.hpp"
9 : #include "pairinteraction/utils/eigen_assertion.hpp"
10 : #include "pairinteraction/utils/eigen_compat.hpp"
11 : #include "pairinteraction/utils/spherical.hpp"
12 :
13 : #include <Eigen/Dense>
14 : #include <algorithm>
15 : #include <limits>
16 : #include <memory>
17 : #include <set>
18 : #include <spdlog/spdlog.h>
19 :
20 : namespace pairinteraction {
21 : template <typename Scalar>
22 134 : SystemAtom<Scalar>::SystemAtom(std::shared_ptr<const basis_t> basis)
23 134 : : System<SystemAtom<Scalar>>(std::move(basis)) {}
24 :
25 : template <typename Scalar>
26 64 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_electric_field(const std::array<real_t, 3> &field) {
27 64 : this->hamiltonian_requires_construction = true;
28 :
29 40 : if (!traits::NumTraits<Scalar>::is_complex_v && field[1] != 0) {
30 0 : throw std::invalid_argument(
31 : "The field must not have a y-component if the scalar type is real.");
32 : }
33 :
34 64 : electric_field = field;
35 :
36 64 : return *this;
37 : }
38 :
39 : template <typename Scalar>
40 24 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_magnetic_field(const std::array<real_t, 3> &field) {
41 24 : this->hamiltonian_requires_construction = true;
42 :
43 20 : if (!traits::NumTraits<Scalar>::is_complex_v && field[1] != 0) {
44 0 : throw std::invalid_argument(
45 : "The field must not have a y-component if the scalar type is real.");
46 : }
47 :
48 24 : magnetic_field = field;
49 :
50 24 : return *this;
51 : }
52 :
53 : template <typename Scalar>
54 15 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_diamagnetism_enabled(bool enable) {
55 15 : this->hamiltonian_requires_construction = true;
56 15 : diamagnetism_enabled = enable;
57 15 : return *this;
58 : }
59 :
60 : template <typename Scalar>
61 : SystemAtom<Scalar> &
62 19 : SystemAtom<Scalar>::set_ion_distance_vector(const std::array<real_t, 3> &vector) {
63 19 : this->hamiltonian_requires_construction = true;
64 :
65 2 : if (!traits::NumTraits<Scalar>::is_complex_v && vector[1] != 0) {
66 0 : throw std::invalid_argument(
67 : "The distance vector must not have a y-component if the scalar type is real.");
68 : }
69 :
70 19 : ion_distance_vector = vector;
71 :
72 19 : return *this;
73 : }
74 :
75 : template <typename Scalar>
76 17 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_ion_charge(real_t charge) {
77 17 : this->hamiltonian_requires_construction = true;
78 17 : ion_charge = charge;
79 17 : return *this;
80 : }
81 :
82 : template <typename Scalar>
83 19 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_ion_interaction_order(int value) {
84 19 : this->hamiltonian_requires_construction = true;
85 19 : if (value < 2 || value > 3) {
86 0 : throw std::invalid_argument("The order of the Rydberg-ion interaction must be 2 or 3");
87 : }
88 19 : ion_interaction_order = value;
89 19 : return *this;
90 : }
91 :
92 : template <typename Scalar>
93 188 : void SystemAtom<Scalar>::construct_hamiltonian() const {
94 188 : auto basis = this->hamiltonian->get_basis();
95 :
96 : // Construct the unperturbed Hamiltonian
97 188 : this->hamiltonian = std::make_unique<OperatorAtom<Scalar>>(basis, OperatorType::ENERGY);
98 188 : this->hamiltonian_is_diagonal = true;
99 188 : bool sort_by_quantum_number_f = true;
100 188 : bool sort_by_quantum_number_m = true;
101 188 : bool sort_by_parity = true;
102 :
103 : // Estimate the numerical precision so that we can decide which terms to keep
104 188 : Eigen::VectorX<real_t> diag = this->hamiltonian->get_matrix().diagonal().real();
105 188 : real_t scale = (diag - diag.mean() * Eigen::VectorX<real_t>::Ones(diag.size())).norm();
106 188 : real_t numerical_precision = 100 * scale * std::numeric_limits<real_t>::epsilon();
107 :
108 188 : real_t typical_magnetic_dipole = 1e2; // ~n^1
109 188 : real_t typical_electric_dipole = 1e4; // ~n^2
110 188 : 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 376 : Eigen::Map<const Eigen::Vector3<real_t>> vector_map(ion_distance_vector.data(),
116 188 : ion_distance_vector.size());
117 188 : real_t distance = vector_map.norm();
118 376 : SPDLOG_DEBUG("Distance to the ion: {}", distance);
119 :
120 : // Dipole order
121 188 : 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 19 : Eigen::Vector3<Scalar> vector_dipole_order =
125 19 : spherical::get_transformator<Scalar>(1) * vector_map / distance;
126 19 : vector_dipole_order = vector_dipole_order.conjugate() / std::pow(distance, 2);
127 :
128 76 : for (int q = -1; q <= 1; ++q) {
129 57 : if (std::abs(vector_dipole_order[q + 1]) * typical_electric_dipole >
130 : numerical_precision) {
131 30 : *this->hamiltonian -= ion_charge * vector_dipole_order[q + 1] *
132 60 : OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_DIPOLE, q);
133 30 : this->hamiltonian_is_diagonal = false;
134 30 : sort_by_quantum_number_f = false;
135 30 : sort_by_parity = false;
136 30 : 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 188 : 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 18 : Eigen::Vector<Scalar, 6> vector_quadrupole_order = spherical::get_transformator<Scalar>(2) *
148 18 : Eigen::KroneckerProduct(vector_map / distance, vector_map / distance);
149 18 : vector_quadrupole_order = 3 * vector_quadrupole_order.conjugate() / std::pow(distance, 3);
150 :
151 108 : for (int q = -2; q <= 2; ++q) {
152 90 : if (std::abs(vector_quadrupole_order[q + 2]) * typical_electric_quadrupole >
153 : numerical_precision) {
154 40 : *this->hamiltonian -= ion_charge * vector_quadrupole_order[q + 2] *
155 80 : OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, q);
156 40 : this->hamiltonian_is_diagonal = false;
157 40 : sort_by_quantum_number_f = false;
158 40 : 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 188 : Eigen::Vector3<Scalar> electric_field_spherical = spherical::get_transformator<Scalar>(1) *
167 188 : Eigen::Map<const Eigen::Vector3<real_t>>(electric_field.data(), electric_field.size());
168 188 : electric_field_spherical = electric_field_spherical.conjugate();
169 :
170 : // Transform the magnetic field to spherical coordinates and take the conjugate
171 188 : Eigen::Vector3<Scalar> magnetic_field_spherical = spherical::get_transformator<Scalar>(1) *
172 188 : Eigen::Map<const Eigen::Vector3<real_t>>(magnetic_field.data(), magnetic_field.size());
173 188 : 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 752 : for (int q = -1; q <= 1; ++q) {
180 564 : if (std::abs(electric_field_spherical[q + 1]) * typical_electric_dipole >
181 : numerical_precision) {
182 119 : *this->hamiltonian -= electric_field_spherical[q + 1] *
183 238 : OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_DIPOLE, q);
184 119 : this->hamiltonian_is_diagonal = false;
185 119 : sort_by_quantum_number_f = false;
186 119 : sort_by_parity = false;
187 119 : 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 752 : for (int q = -1; q <= 1; ++q) {
196 564 : if (std::abs(magnetic_field_spherical[q + 1]) * typical_magnetic_dipole >
197 : numerical_precision) {
198 22 : *this->hamiltonian -= magnetic_field_spherical[q + 1] *
199 44 : OperatorAtom<Scalar>(basis, OperatorType::MAGNETIC_DIPOLE, q);
200 22 : this->hamiltonian_is_diagonal = false;
201 22 : sort_by_quantum_number_f = false;
202 22 : 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 188 : if (diamagnetism_enabled) {
214 7 : if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
215 : numerical_precision) {
216 5 : *this->hamiltonian += static_cast<real_t>(1 / 12.) *
217 5 : static_cast<Scalar>(std::pow(magnetic_field_spherical[1], 2)) *
218 10 : OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, 0);
219 5 : *this->hamiltonian -= static_cast<real_t>(1 / 12.) *
220 5 : static_cast<Scalar>(std::pow(magnetic_field_spherical[1], 2)) *
221 10 : OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, 0);
222 5 : this->hamiltonian_is_diagonal = false;
223 5 : sort_by_quantum_number_f = false;
224 : }
225 7 : if (std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
226 9 : numerical_precision &&
227 2 : std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
228 : numerical_precision) {
229 2 : *this->hamiltonian -= static_cast<real_t>(2 / 12.) * magnetic_field_spherical[0] *
230 2 : magnetic_field_spherical[2] *
231 4 : OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, 0);
232 2 : *this->hamiltonian -= static_cast<real_t>(1 / 12.) * magnetic_field_spherical[0] *
233 2 : magnetic_field_spherical[2] *
234 4 : OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, 0);
235 2 : this->hamiltonian_is_diagonal = false;
236 2 : sort_by_quantum_number_f = false;
237 : }
238 7 : if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
239 12 : numerical_precision &&
240 5 : std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
241 : numerical_precision) {
242 0 : *this->hamiltonian -= static_cast<real_t>(std::sqrt(3.0) / 12.) *
243 0 : magnetic_field_spherical[1] * magnetic_field_spherical[2] *
244 0 : OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, 1);
245 0 : this->hamiltonian_is_diagonal = false;
246 0 : sort_by_quantum_number_f = false;
247 0 : sort_by_quantum_number_m = false;
248 : }
249 7 : if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
250 12 : numerical_precision &&
251 5 : std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
252 : numerical_precision) {
253 0 : *this->hamiltonian -= static_cast<real_t>(std::sqrt(3.0) / 12.) *
254 0 : magnetic_field_spherical[1] * magnetic_field_spherical[0] *
255 0 : OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, -1);
256 0 : this->hamiltonian_is_diagonal = false;
257 0 : sort_by_quantum_number_f = false;
258 0 : sort_by_quantum_number_m = false;
259 : }
260 7 : if (std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
261 : numerical_precision) {
262 2 : *this->hamiltonian -= static_cast<real_t>(std::sqrt(1.5) / 12.) *
263 2 : static_cast<Scalar>(std::pow(magnetic_field_spherical[2], 2)) *
264 4 : OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, 2);
265 2 : this->hamiltonian_is_diagonal = false;
266 2 : sort_by_quantum_number_f = false;
267 2 : sort_by_quantum_number_m = false;
268 : }
269 7 : if (std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
270 : numerical_precision) {
271 2 : *this->hamiltonian -= static_cast<real_t>(std::sqrt(1.5) / 12.) *
272 2 : static_cast<Scalar>(std::pow(magnetic_field_spherical[0], 2)) *
273 4 : OperatorAtom<Scalar>(basis, OperatorType::ELECTRIC_QUADRUPOLE, -2);
274 2 : this->hamiltonian_is_diagonal = false;
275 2 : sort_by_quantum_number_f = false;
276 2 : sort_by_quantum_number_m = false;
277 : }
278 : }
279 :
280 : // Store which labels can be used to block-diagonalize the Hamiltonian
281 188 : this->blockdiagonalizing_labels.clear();
282 188 : if (sort_by_quantum_number_f) {
283 95 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_F);
284 : }
285 188 : if (sort_by_quantum_number_m) {
286 143 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_M);
287 : }
288 188 : if (sort_by_parity) {
289 114 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_PARITY);
290 : }
291 188 : }
292 :
293 : // Explicit instantiations
294 : template class SystemAtom<double>;
295 : template class SystemAtom<std::complex<double>>;
296 : } // namespace pairinteraction
|
