diff --git a/examples/fei2_tutorial.jl b/examples/fei2_tutorial.jl index 616e7d3fb..a8a83bbf7 100644 --- a/examples/fei2_tutorial.jl +++ b/examples/fei2_tutorial.jl @@ -211,7 +211,7 @@ print_wrapped_intensities(sys) # ground state. # If the desired ground state is already known, as with FeI$_2$, it could be -# entered by hand using [`polarize_spin!`](@ref). Alternatively, in the case of +# entered by hand using [`set_dipole!`](@ref). Alternatively, in the case of # FeI$_2$, we could repeatedly employ the above randomization and minimization # procedure until a defect-free configuration is found. Some systems will have # more complicated ground states, which can be much more challenging to find. @@ -245,14 +245,12 @@ plot_spins(sys_min; arrowlength=2.5, linewidth=0.75, arrowsize=1.5) swt = SpinWaveTheory(sys_min); -# The dispersion relation can be calculated by providing `dispersion` with a -# `SpinWaveTheory` and a list of wave vectors. For each wave vector, -# `dispersion` will return a list of energies, one for each band. -# Frequently one wants dispersion information along lines -# that connect special wave vectors. The function [`connected_path`](@ref) -# linearly samples between provided $q$-points, with a given sample density. - -points = [[0, 0, 0], # List of wave vectors that define a path +# Next, select a sequence of wavevectors that will define a piecewise linear +# path in reciprocal space. The coefficients of each $q$-vector are given in +# reciprocal lattice units (RLU), i.e., are multiples of the recripocal lattice +# vectors. The function [`connected_path`](@ref) will linearly sample between +# the provided $q$-points with some `density`. +points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [1/2, 0, 0], @@ -262,21 +260,32 @@ labels = string.(points) density = 600 path, markers = connected_path(swt, points, density); -# `dispersion` may now be called on the wave vectors along the generated path. -# Each column of the returned matrix corresponds to a different mode. +# The dispersion relation can be calculated by providing `dispersion` with a +# `SpinWaveTheory` and a list of wave vectors. For each wave vector, +# `dispersion` will return a list of energies, one for each band. Frequently one +# wants dispersion information along lines that connect special wave vectors. +# The function [`connected_path`](@ref) linearly samples between provided +# $𝐪$-points, with a given sample density. + +# The [`dispersion`](@ref) function defines the quasiparticle excitation +# energies $ω_i(𝐪)$ for each point $𝐪$ along the reciprocal space path. disp = dispersion(swt, path); -# In addition to the band energies ``\omega_i``, Sunny can calculate -# the inelastic neutron scattering intensity ``I(q,\omega_i(q))`` according to -# an [`intensity_formula`](@ref). The default formula applies a polarization correction ``(1 - Q\otimes Q)``. +# In addition to the band energies $ω_i(𝐪)$, Sunny can calculate the inelastic +# neutron scattering intensity $I_i(𝐪)$ for each band $i$ according to an +# [`intensity_formula`](@ref). The default formula applies a polarization +# correction $(1 - 𝐪⊗𝐪)$. Selecting `delta_function_kernel` specifies that we +# want the energy and intensity of each band individually. + +formula = intensity_formula(swt; kernel=delta_function_kernel) -formula = intensity_formula(swt; kernel = delta_function_kernel) +# The function [`intensities_bands`](@ref) uses linear spin wave theory to +# calculate both the dispersion and intensity data for the provided path. -# The `delta_function_kernel` specifies that we want the energy and intensity of each -# band individually. +disp, intensity = intensities_bands(swt, path; formula); -disp, intensity = intensities_bands(swt, path; formula) +# These can be plotted in GLMakie. fig = Figure() ax = Axis(fig[1,1]; xlabel="𝐪", ylabel="Energy (meV)", @@ -288,44 +297,78 @@ for i in axes(disp)[2] end fig -# For comparison with inelastic neutron scattering (INS) data, it's often -# appropriate to apply a lorentzian broadening to the bands. +# To make comparisons with inelastic neutron scattering (INS) data, it is helpful to employ +# an empirical broadening kernel. + +γ = 0.15 # width in meV +broadened_formula = intensity_formula(swt; kernel=lorentzian(γ)) + +# The broadened intensity can be calculated for using +# [`intensities_broadened`](@ref). One must space a path in $𝐪$-space and an +# energy range. -γ = 0.05 # Lorentzian broadening parameter -broadened_formula = intensity_formula(swt; kernel = lorentzian(γ)) +energies = collect(0:0.01:7.5) # 0 < ω < 7.5 (meV). +is1 = intensities_broadened(swt, path, energies, broadened_formula); -# The broadened intensity can be calculated for any ``(q,\omega)`` -# using `intensities_broadened`: -energies = collect(0.0:0.01:7.5) # Energies to calculate -is = intensities_broadened(swt, path, energies, broadened_formula) +# In real FeI$_2$, there will be competing magnetic domains associated with the +# three possible orientations of the symmetry-broken ground state, which derive +# from the three-fold rotational symmetry of a triangular lattice. +# +# The matrix $A$ defined below implements a 120° clockwise rotation on +# wavevectors $𝐪$ defined in RLU. + +A = [-1 -1 0; + 1 0 0; + 0 0 1] +@assert A^3 ≈ [1 0 0; 0 1 0; 0 0 1] # Rotation by 360° is the identity + +# We will average the intensity data over the three symmetry-equivalent +# rotations of the $𝐪$-space path. + +function rotated_intensity(R) + rotated_path = [R * q for q in path] + intensities_broadened(swt, rotated_path, energies, broadened_formula) +end + +@assert is1 == rotated_intensity(A^0) +is2 = rotated_intensity(A) +is3 = rotated_intensity(A^2) +is_averaged = (is1 + is2 + is3) / 3 fig = Figure() ax = Axis(fig[1,1], xlabel="(H,0,0)", ylabel="Energy (meV)", xticks=(markers, labels), - xticklabelrotation=π/8) -heatmap!(ax, 1:size(is, 1), energies, is) + ) +heatmap!(ax, 1:size(is_averaged, 1), energies, is_averaged) fig -# The existence of a lower-energy, single-ion bound state is in qualitative -# agreement with the experimental data in [Bai et -# al.](https://doi.org/10.1038/s41567-020-01110-1) (Note that the publication -# figure uses a different coordinate system to label the same wave vectors and -# the experimental data necessarily averages over the three degenerate ground -# states.) +# This result can be directly compared to experimental neutron scattering data +# [Bai et al.](https://doi.org/10.1038/s41567-020-01110-1) # ```@raw html # # ``` +# +# (Note that the publication figure uses a different coordinate system to label +# the wave vectors.) +# +# To get this agreement, the use of SU(3) coherent states is essential (in other +# words, we needed a theory of multi-flavored bosons). The lower band has large +# quadrupolar character, and arises from the strong easy-axis anisotropy of +# FeI$_2$. By setting `mode = :SUN`, the calculation captures this coupled +# dipole-quadrupole dynamics. +# +# An interesting exercise is to repeat the same study, but using `mode = +# :dipole` instead of `:SUN`. That alternative choice would constrain the +# coherent state dynamics to the space of dipoles only. -# The full data from the dynamical spin structure factor (DSSF) can be -# retrieved with the `dssf` function. Like `dispersion` and `intensities`, -# `dssf` takes an array of wave vectors. +# The full dynamical spin structure factor (DSSF) can be retrieved as a $3×3$ +# matrix with the [`dssf`](@ref) function, for a given path of $𝐪$-vectors. -disp, Sαβ = dssf(swt, [[0, 0, 0]]); +disp, is = dssf(swt, path); -# `disp` is identical to the output that is obtained from `dispersion` and -# contains the energy of each mode at the specified wave vectors. `Sαβ` contains -# a $3×3$ matrix for each of these modes. For example, `Sαβ[1]` are the -# dynamical structure factor intensities for the first mode, and `Sαβ[1][2,3]` -# is associated with the $α=ŷ$ and $β=ẑ$ spin components. +# The first output `disp` is identical to that obtained from `dispersion`. The +# second output `is` contains a list of $3×3$ matrix of intensities. For +# example, `is[q,n][2,3]` yields the $(ŷ,ẑ)$ component of the structure factor +# intensity for `nth` mode at the `q`th wavevector in the `path`. # ## What's next? #