phase diagram of 3D Anderson model

L = 10; % system size

nr = 10;

nsteps = 100;

W_values = linspace(0, 20, nsteps);

IPR_over_W = cell(1, nsteps);

eigenvalues_cell = cell(nr, nsteps);

for W_index = 1:nsteps

W = W_values(W_index);

IPRs_for_this_W = zeros(1, nr);

for ir = 1:nr

H = generate_hamiltonian(L, W);

[eigenstates, D] = eig(H);

eigenvalues = diag(D);

energy_levels{ir} = eigenvalues(:)';

IPR{ir} = compute_IPR(eigenstates);

IPRs_for_this_W(ir) = mean(IPR{ir});

% Store the eigenvalues in a cell array

eigenvalues_cell{ir, W_index} = eigenvalues(:)';

end

% Store the mean IPR over 100 realizations for this W value

IPR_over_W{W_index} = mean(IPRs_for_this_W);

disp([W, IPR_over_W{W_index}])

end

%

% Create a figure for the color plot

figure;

% Define the colormap and use the IPR values for coloring

cmap = jet;

% Loop over W values

for W_index = 1:nsteps

W = W_values(W_index);

% Extract eigenvalues from the cell array

eigenvalues_for_W = eigenvalues_cell(:, W_index);

% Initialize arrays to store eigenvalues and y-axis values

eigenvalues_numeric = [];

y_values = [];

cdata = []; % Initialize color data

% Convert eigenvalues_for_W to a numeric array and create y_values

for ir = 1:nr

eigenvalues_ir = eigenvalues_for_W{ir};

eigenvalues_numeric = [eigenvalues_numeric, eigenvalues_ir];

y_values = [y_values, W * ones(size(eigenvalues_ir))];

% Calculate IPR values (You should have a function for this)

IPR_value = compute_IPR(eigenvalues_ir);

cdata = [cdata, IPR_value];

end

% Create a scatter plot with eigenvalues on the x-axis and y_values on the y-axis

scatter(eigenvalues_numeric, y_values, 20, cdata*1e-4, 'filled');

% imagesc(W_values, 1:nr, IPR_matrix);

colormap(cmap);

colorbar; % Add a colorbar

hold on;

end

xlabel('Eigenvalues');

ylabel('W');

title('Color Plot: Eigenvalues vs. W vs. IPR');

hold off;

%%% box disorder

function disorder = generate_disorder(L, W)

disorder = W * ((rand(L, L, L) - 0.5));

end

% Generate the Hamiltonian matrix for one realization of the random disorder

function H = generate_hamiltonian(L, W)

H = zeros(L * L * L, L * L * L);

disorder = generate_disorder(L, W);

for i = 1:L

ip1 = mod(i, L) + 1;

for j = 1:L

jp1 = mod(j, L) + 1;

for k = 1:L

kp1 = mod(k, L) + 1;

H((i - 1) * L * L + (j - 1) * L + k, (i - 1) * L * L + (j - 1) * L + k) = disorder(i, j, k);

H(ip1 * L * L + j * L + k, (i - 1) * L * L + j * L + k) = 1.0;

H((i - 1) * L * L + jp1 * L + k, (i - 1) * L * L + j * L + k) = 1.0;

H((i - 1) * L * L + j * L + kp1, (i - 1) * L * L + j * L + k) = 1.0;

H((i - 1) * L * L + j * L + k, ip1 * L * L + j * L + k) = 1.0;

H((i - 1) * L * L + j * L + k, (i - 1) * L * L + jp1 * L + k) = 1.0;

H((i - 1) * L * L + j * L + k, (i - 1) * L * L + j * L + kp1) = 1.0;

end

end

end

end

function IPR = compute_IPR(eigenstates)

IPR = sum(abs(eigenstates).^4, 1);

end

Answer 1

[Abhinaya Kennedy]

Hi Zahid,

As per my understanding, you are trying to visualize the transition from extended to localized states in the Anderson model by plotting the inverse participation ratio (IPR) across different disorder strengths (W) and eigenvalues.

The issue you're facing with the center of the band not showing the expected transition could be due to several reasons:

1. Color Mapping Range: The colormap might not be reflecting the changes in IPR values because the range of the colormap is not set correctly. If the color range is too broad, small changes in IPR might not be visible.

2. IPR Calculation: There might be an issue with how the IPR is being calculated or averaged. Ensure that the `compute_IPR` function is correctly implemented and is accurately reflecting the localization of states.

3. Numerical Resolution: The discretization of W values and the number of realizations might not be enough to capture the transition accurately. Increasing the `nsteps` and `nr` could help.

4. Visualization Technique: The method of visualization might not be sensitive enough to show the subtle changes in IPR. Consider using a different visualization approach or adjusting the parameters of the scatter plot.

Looking at the code, I see a potential issue in the loop where you plot the scatter points. You calculate `IPR_value` using `compute_IPR(eigenvalues_ir)`, but it seems like you should be using the already computed IPR values which are stored in `IPR{ir}` instead of recalculating them. This could be a mistake since `eigenvalues_ir` is just a list of eigenvalues, not the eigenstates required for the IPR calculation.

Here's the corrected section of the code:

% ...
% Loop over W values
for W_index = 1:nsteps
    W = W_values(W_index);
    % Extract eigenvalues from the cell array
    eigenvalues_for_W = eigenvalues_cell(:, W_index);
    % Extract IPR from the cell array
    IPR_for_W = IPR_over_W{W_index};
    % Initialize arrays to store eigenvalues and y-axis values
    eigenvalues_numeric = [];
    y_values = [];
    cdata = []; % Initialize color data
    % Convert eigenvalues_for_W to a numeric array and create y_values
    for ir = 1:nr
        eigenvalues_ir = eigenvalues_for_W{ir};
        eigenvalues_numeric = [eigenvalues_numeric, eigenvalues_ir];
        y_values = [y_values, W * ones(size(eigenvalues_ir))];
        % Use the precomputed IPR values
        cdata = [cdata, IPR{ir}];
    end
    % Create a scatter plot with eigenvalues on the x-axis and y_values on the y-axis
    scatter(eigenvalues_numeric, y_values, 20, cdata, 'filled'); % Removed the scaling of IPR values by 1e-4
    colormap(cmap);
    colorbar; % Add a colorbar
    hold on;
end
% ...

Make sure that the `cdata` array contains the correct IPR values for each eigenvalue and that the color scaling in the scatter plot reflects the range of IPR values you expect to see.

Additionally, consider setting the limits of the color axis (`caxis`) to a range that emphasizes the changes in IPR values around the critical disorder strength. This will help you visualize the transition more clearly.

If you continue to face issues, it would be helpful to plot the raw IPR data as a function of W for a few eigenvalues to manually inspect if the transition is occurring as expected. If the raw data shows the transition, then the issue is likely with the visualization and not the computation.

Hope this helps!