# Numerical Solution to Schrödinger Equation--Multiple Wells

I am trying to solve for the allowed wavefunctions and energies for a 1D quartic potential well.

To do this I am using the patching method (https://engineering.dartmouth.edu/microeng/otherweb/henning/papers/jap-sub95.pdf). I integrate from x = 0 to x = patch and from x = patch to x = end with the Numerov method. I then find the zeros of the function

$$\psi_{left}(patch) *\frac{\mathrm{d} \psi_{right}(patch) }{\mathrm{d} x} - \psi_{right}(patch) *\frac{\mathrm{d} \psi_{left}(patch) }{\mathrm{d} x}$$

where left means integration from x = 0 to x = patch, and right means integration from x = end to x = patch.

I am able to find the correct wavefunctions for a quartic well of small depth. The potential: And the first four wavefunctions that were found: The issue I am having is that when I increase the well depth, I am unable to find the correct wavefunctions. For the following deeper well: The first found wavefunction (below) for this deeper well is clearly not the first allowed energy because it has too many zero crossings: When I inspect the minimizing function that I find the roots of, I see that it does not cross zero at low energies.

I suspect the issue for this is that in the deeper well, the two potential wells can be considered separately because there is so little wavefunction overlap. I would, however, like my solver to be robust enough to find the isolated wavefunctions, or at the very least recognize when the two potential wells can be considered separately. I would also like my solver to solve asymmetric wells, which I am not sure it can do right now. Is there a way I can adapt my code to solve these two issues?

Also, here is the MATLAB code reproduced below. The main function is the first code block and is the script.

%% Outputs
% E_allowed; % Allowed eigenvalues
% wavefunctions_allowed; % Allowed eigenfunctions

%% Inputs
potential_function; % potential function, energy in eV
z_axis; % the 1D problem
h; % z_step, z(2) - z(1), distance in meters
mass_particle;
E_guess_max; E_guess_min;
E_guess_step;
E_disp;

%% Set up
% Find eigenfunctions and eigenenergies for an arbitrary potential function
if(~any(strcmp(who,'MAX_SHRO_LOOP')))
MAX_SHRO_LOOP = 200; % Maximum runs for the bisection loop
end

if(~any(strcmp(who,'TWO_PASS_PATCHING')))
TWO_PASS_PATCHING = 1; % Maximum runs for the bisection loop
end

if(~any(strcmp(who,'E_guess_step')))
E_guess_step = 1e-4; % MUST BE SET ACCORDINGLY
%MUST BE SMALLER THAN THE SPACING BETWEEN ADJACENET ENERGY LEVELS
%MUST BE LARGE SUCH THAT WE INCREMENT E AFTER EIGEN FOUND WE SEE CHANGE
end

if(~any(strcmp(who,'z_patch_vec')))
% the x patches that can be tried
z_patch_vec = 3 : numel(z_axis) - 2;
end

if(~any(strcmp(who,'z_patch_index')))
z_patch_index = round(numel(z_patch_vec) / 2); % two pass z patching
end

if(~any(strcmp(who,'init_small')))
init_small = 1e-100; % have to adjust to get good results
% if determinant is too small then increase this number--if the eigenvalues
% are found too quickly
% if NaNs in wavefunctions then decrease
end

if(~any(strcmp(who,'num_to_find')))
num_to_find = [30 60 100 500 Inf];
end

if(~any(strcmp(who,'accuracy_to_find')))
accuracy_to_find = [1e-16 1e-8 1e-8 1e-8 1e-8];
% after each num_to_find has been reached, change the accuracy
end

% E in eV
if(~any(strcmp(who,'E_guess_min')))
E_guess_min = min(potential_function);
end

if(~any(strcmp(who,'E_guess_max')))
E_guess_max = max([potential_function(1) potential_function(end)]);
end

if(~any(strcmp(who,'E_disp')))
E_disp = 0;
end

%%
mass_particle = mass_particle / q;

Es_guessed = [];
crossfires_guessed = [];

% the calculated stuff...
E_allowed = [];
wavefunctions_allowed = [];
left_allowed_wavefunctions = [];
right_allowed_wavefunctions = [];

E_guess = E_guess_min;
E_bracket_low = E_guess_min;
E_bracket_high = E_guess_min;

% for displaying stuff
E_guess_range = abs(E_guess_max - E_guess_min);
curr_percent = 10;
fprintf('\nShrodinger: ');

z_patches = []; % the x_patches tried
if_loop_counter = 1; % the counter for the adaptive x pacthing

while E_guess < E_guess_max

EIGEN_ERR = accuracy_to_find(sum(numel(E_allowed) > num_to_find) + 1);
z_patch = z_patch_vec(z_patch_index);
% print progress bar
if( 100 * abs(min(potential_function) - E_guess) / E_guess_range >= curr_percent )
fprintf(1,'%2d%%... ',(curr_percent));
curr_percent = curr_percent + 10;
end

% Find a good x_patch value

% Low energy bracket
T = h^2 * 2 * mass_particle * (potential_function - E_bracket_low) / (12 * h_bar^2); % T function
guess_left_bracket_low  = Numerov_Left(0, init_small, T);
guess_right_bracket_low = Numerov_Right(0, init_small, T);
[crossfire_mat_bracket_low val_low] = make_crossfire_mat(guess_left_bracket_low, guess_right_bracket_low, z_patch, h);

% high energy bracket
E_bracket_high = E_bracket_low + E_guess_step;
T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function
guess_left_bracket_high  = Numerov_Left(0, init_small, T);
guess_right_bracket_high = Numerov_Right(0, init_small, T);
crossfire_mat_bracket_high = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);
%--this code can be deleted...

% poor patching crossfire returns NaN, try another patching value
if(isnan(crossfire_mat_bracket_low) ||  isnan(crossfire_mat_bracket_high))

if_loop_counter = 1 + if_loop_counter;
z_patch_index = mod(z_patch_index + 1, numel(z_patch_vec)) + 1;
z_patch = z_patch_vec(z_patch_index);

if(if_loop_counter > numel(z_patch_vec))
disp('patching'); % all the x patchings were tried, none worked...
break;
end

end
%--this code can be deleted...

%--this if statement (but not the code in the block) can be deleted
% only continue if the x patching is okay and neither crossfire is not NaN
if(~isnan(crossfire_mat_bracket_low) && ~isnan(crossfire_mat_bracket_high))

if_loop_counter = 1; % reset the x patch counter so we can do it again later if needed
z_patches = [z_patches; z_patch]; % the x_patches we have tried

crossfires_guessed = [crossfires_guessed; val_low];
Es_guessed = [Es_guessed; E_bracket_low];

%-------------------------- bracket guesses --------------------------%
% Want to bound the eigenenerngy between the low and high brackets
% Now having found a good x_patch, let us continue to bracket the eigenvalue guess

% low energy bracket is either the min energy (if 1st try) or the
% previous eigenenergy + guess_step

% high energy bracket search
% determinant has to flip signs
while(crossfire_mat_bracket_low * crossfire_mat_bracket_high > 0)

E_bracket_high = E_bracket_high + E_guess_step;

T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function

guess_left_bracket_high  = Numerov_Left(0, init_small, T);
guess_right_bracket_high = Numerov_Right(0, init_small, T);

[crossfire_mat_bracket_high val] = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);
crossfires_guessed = [crossfires_guessed; val];
Es_guessed = [Es_guessed; E_bracket_high];

end

counter = 1; % counter for the binary search
found_flag = 0; % 1 = found an eigenenergy

% set the first guess to be the low bracket
E_guess_previous = E_bracket_low;
crossfire_mat_previous = crossfire_mat_bracket_low;

% guess
E_guess = (E_bracket_high + E_bracket_low) / 2;

% binary search
% search until MATLAB resolution reached or iteration counter satisfied
while(abs(E_guess_previous - E_guess) > EIGEN_ERR && counter < MAX_SHRO_LOOP && ~found_flag)

if(E_guess > E_guess_max)
found_flag = 1;
break
end

% for the bisection guess
T = h^2 * 2 * mass_particle * (potential_function - E_guess) / (12 * h_bar^2); % T function
guess_left  = Numerov_Left(0, init_small, T);
guess_right = Numerov_Right(0, init_small, T);
[crossfire_mat val] = make_crossfire_mat(guess_left, guess_right, z_patch, h);

crossfires_guessed = [crossfires_guessed; val];
Es_guessed = [Es_guessed; E_guess];

% E guess correct
if(crossfire_mat == 0)

E_allowed = [E_allowed; E_guess];
right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];
left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];
found_flag = 1; % will not add to the allowed stuff at the end of the bisection loop

% guessed too low, no sign change
elseif(crossfire_mat * crossfire_mat_previous > 0)

E_bracket_low = E_guess;

else % guessed too high and the sign changed

E_bracket_high = E_guess;

end

% diagnostic display
if(E_disp)
disp(E_guess);
end

% go onto the next bisection loop iteration
E_guess_previous = E_guess;
crossfire_mat_previous = crossfire_mat;
E_guess = 0.5 * (E_bracket_low + E_bracket_high); % New guess

counter = counter + 1;

end

% if the binary search stopped
% stopped because the difference between two eigenenergies was so
% small or because the max iterations was reached...not because
% determinant was zero
if(~found_flag)

E_allowed = [E_allowed; E_guess];
right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];
left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];

end

% increment the brackets for the next run
E_bracket_low = E_allowed(end) + E_guess_step;
E_bracket_high = E_bracket_low + E_guess_step;

end

end

%%
% Clean up wavefunctions
% wavefunction correction
wavefunctions_allowed = zeros(numel(z_axis), numel(E_allowed));
patching = [];
for i = 1 : numel(E_allowed)

[wavefunctions_allowed(:,i) patch_out] = wavefunction_patch(left_allowed_wavefunctions(:,i), right_allowed_wavefunctions(:,i), z_axis, z_patches(i), E_allowed(i), potential_function);
patching = [patching; patch_out];

end

% flag_zerocross_eigenlevel = 1; % so far ok
% % num nodes ~= the energy level
% for i = 1 : numel(E_allowed)
%
%     % zero crossings not commensurate with eigenvalue level
%     if(numel(find_roots(wavefunctions_allowed(:,i))) ~= i - 1)
%         flag_zerocross_eigenlevel = 0;
%         break
%     end
%
% end
%
% if(~flag_zerocross_eigenlevel)
%
%     %disp('ERROR');
%     % the wells are too seperate, consider them as isolated!!
%
% end

%% Two pass patching, now we do it again...having found the optimal patching where the agreement is closest. Necessary?
% the calculated stuff...
if(TWO_PASS_PATCHING)

E_allowed = [];
wavefunctions_allowed = [];
left_allowed_wavefunctions = [];
right_allowed_wavefunctions = [];

% E in eV
E_guess_min = min(potential_function);
E_guess_max = max([potential_function(1) potential_function(end)]);

E_guess = E_guess_min;
E_bracket_low = E_guess_min;
E_bracket_high = E_guess_min;

% for displaying stuff
E_guess_range = abs(E_guess_max - E_guess_min);
curr_percent = 10;
fprintf('\nShrodinger: ');

while E_guess < E_guess_max && numel(E_allowed) < numel(patching)

EIGEN_ERR = accuracy_to_find(sum(numel(E_allowed) > num_to_find) + 1);

% second pass patching
z_patch = patching(numel(E_allowed) + 1);

% print progress bar
if( 100 * abs(min(potential_function) - E_guess) / E_guess_range >= curr_percent )
fprintf(1,'%2d%%... ',(curr_percent));
curr_percent = curr_percent + 10;
end

% Low energy bracket
T = h^2 * 2 * mass_particle * (potential_function - E_bracket_low) / (12 * h_bar^2); % T function
guess_left_bracket_low  = Numerov_Left(0, init_small, T);
guess_right_bracket_low = Numerov_Right(0, init_small, T);
[crossfire_mat_bracket_low val_low] = make_crossfire_mat(guess_left_bracket_low, guess_right_bracket_low, z_patch, h);

% high energy bracket
E_bracket_high = E_bracket_low + E_guess_step;
T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function
guess_left_bracket_high  = Numerov_Left(0, init_small, T);
guess_right_bracket_high = Numerov_Right(0, init_small, T);
crossfire_mat_bracket_high = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);

crossfires_guessed = [crossfires_guessed; val_low];
Es_guessed = [Es_guessed; E_bracket_low];

%-------------------------- bracket guesses --------------------------%
% Want to bound the eigenenerngy between the low and high brackets
% Now having found a good x_patch, let us continue to bracket the eigenvalue guess

% low energy bracket is either the min energy (if 1st try) or the
% previous eigenenergy + guess_step

% high energy bracket search
% determinant has to flip signs
while(crossfire_mat_bracket_low * crossfire_mat_bracket_high > 0)

E_bracket_high = E_bracket_high + E_guess_step;

T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function

guess_left_bracket_high  = Numerov_Left(0, init_small, T);
guess_right_bracket_high = Numerov_Right(0, init_small, T);

[crossfire_mat_bracket_high val] = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);
crossfires_guessed = [crossfires_guessed; val];
Es_guessed = [Es_guessed; E_bracket_high];

end

counter = 1; % counter for the binary search
found_flag = 0; % 1 = found an eigenenergy

% set the first guess to be the low bracket
E_guess_previous = E_bracket_low;
crossfire_mat_previous = crossfire_mat_bracket_low;

% guess
E_guess = (E_bracket_high + E_bracket_low) / 2;

% binary search
% search until MATLAB resolution reached or iteration counter satisfied
while(abs(E_guess_previous - E_guess) > EIGEN_ERR && counter < MAX_SHRO_LOOP && ~found_flag)

if(E_guess > E_guess_max)
found_flag = 1;
break
end

% for the bisection guess
T = h^2 * 2 * mass_particle * (potential_function - E_guess) / (12 * h_bar^2); % T function
guess_left  = Numerov_Left(0, init_small, T);
guess_right = Numerov_Right(0, init_small, T);
[crossfire_mat val] = make_crossfire_mat(guess_left, guess_right, z_patch, h);

crossfires_guessed = [crossfires_guessed; val];
Es_guessed = [Es_guessed; E_guess];

% E guess correct
if(crossfire_mat == 0)

E_allowed = [E_allowed; E_guess];
right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];
left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];
found_flag = 1; % will not add to the allowed stuff at the end of the bisection loop

% guessed too low, no sign change
elseif(crossfire_mat * crossfire_mat_previous > 0)

E_bracket_low = E_guess;

else % guessed too high and the sign changed

E_bracket_high = E_guess;

end

% diagnostic display
if(E_disp)
disp(E_guess);
end

% go onto the next bisection loop iteration
E_guess_previous = E_guess;
crossfire_mat_previous = crossfire_mat;
E_guess = 0.5 * (E_bracket_low + E_bracket_high); % New guess

counter = counter + 1;

end

% if the binary search stopped
% stopped because the difference between two eigenenergies was so
% small or because the max iterations was reached...not because
% determinant was zero
if(~found_flag)

E_allowed = [E_allowed; E_guess];
right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];
left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];

end

% increment the brackets for the next run
E_bracket_low = E_allowed(end) + E_guess_step;
E_bracket_high = E_bracket_low + E_guess_step;

end

%%
% Clean up wavefunctions
% wavefunction correction
wavefunctions_allowed = zeros(numel(z_axis), numel(E_allowed));
patching_second_pass = [];
for i = 1 : numel(E_allowed)

[wavefunctions_allowed(:,i) patch_out] = wavefunction_patch(left_allowed_wavefunctions(:,i), right_allowed_wavefunctions(:,i), z_axis, patching(i), E_allowed(i), potential_function);
patching_second_pass = [patching_second_pass; patch_out];

end

end

clearvars MAX_SHRO_LOOP TWO_PASS_PATCHING E_guess_step z_patch_vec z_patch_index init_small...
accuracy_to_find num_to_find E_guess_min E_guess_max


function u = Numerov_Left(init1, init2, T)
% init1 = x(1) wavefunction value, must be zero for bound particle
% init2 = x(2) wavefunction value, can be anything b/c normalization
% T function for Numerov

u = zeros(numel(T), 1);
u(1) = init1; % boundary condition
u(2) = init2; % boundary condition

for i = 2 : numel(T) - 1
% Solve for u(x + h)
u(i + 1) = (( (2 + 10 * T(i)) * u(i) - (1 - T(i - 1)) * u(i - 1) ) / (1 - T(i + 1)));

end

end


function [sign_val, val] = make_crossfire_mat(left, right, patch_index, h)

crossfire = [left(patch_index) right(patch_index);...
(left(patch_index) - left(patch_index - 1)) / h ...
(right(patch_index) - right(patch_index -1)) / h];
val = det(crossfire);

sign_val = sign(val);

end


function u = Numerov_Right(init1, init2, T)
% init1 = x(1) wavefunction value, must be zero for bound particle
% init2 = x(2) wavefunction value, can be anything b/c normalization
% T function for Numerov

u = zeros(numel(T), 1);
u(end) = init1; % boundary condition
u(end - 1) = init2; % boundary condition

for i = numel(T) - 1 : -1 : 2
% Solve for u(x - h)
u(i - 1) = (( (2 + 10 * T(i)) * u(i) - (1 - T(i + 1)) * u(i + 1) ) / (1 - T(i - 1)));

end

end


function [out patch] = wavefunction_patch(left, right, z_axis, z_patch, eigenenergy, band_function)
% patch the left and right wavefunctions together

% for symmemtric potentials
if(sum(band_function - fliplr(flipud(band_function))) == 0)
patch = round(numel(band_function) / 2);
else

% patch at location that is closest and nonzero
% the division will make it nonzero because zero will go to inf!!
left = left * right(z_patch) / left(z_patch);
[val patch] = min((abs(left.^2 - right.^2)) ./ (left.^2));

% patch at a peak
patch = temp(1); % return the patching point so we can go back and rerun the crossfire with this new patching point
end

out = [left(1 : patch); right(patch + 1 : end)];
out = out / sqrt(trapz(z_axis, out.^2));

end

• Do you need to use Numerov's method? – nicoguaro Aug 5 '16 at 2:36
• No, do you think an alternate integration scheme would be better? – Geoffrey Xiao Aug 5 '16 at 4:56
• It would be a good idea if you write the potential for your equation and the figures of your eigenvalues. I don't know about this method, that is why I asked. But if you can use other methods like Finite Differences, Finite Elements or Ritz method. I discussed a little bit in this answer. – nicoguaro Aug 8 '16 at 20:42
• These kinds of questions are usually much easier to answer if you post your code directly (if possible), not just a description of the method in a paper (which isn't usually sufficient). Is it possible to see the code? – Kirill Aug 9 '16 at 1:20
• A standard method is to use basis set expansion. Problems of the type you describe can be solved easily by using particle in a box wavefunctions as a basis, see J. Chem. Education. 2011, vol 88, p 929 for examples and sone code. – porphyrin Aug 9 '16 at 19:14

While I cannot help you with your specific implementation, I want to point out to an alternative method (as already indicated in a comment to phil's answer) : Marston's "Fourier Grid Hamiltonian" (FGH).

The FGH is a pseudo-spectral method and yields a simple recipe for constructing the discretized Hamiltonian matrix for bound systems. As usual, eigenvectors and eigenvalues of the discrete Hamiltonian represent discrete wavefunctions and corresponding energies. For the FGH, the Hamiltonian $H = T+V$ is split in a kinetic energy part $T$ and the potential part $V$. In position space the matrix elements $\left< x\vert V\vert x'\right>$ are just $V(x)\delta(x- x')$, so the potential energy contributes only the Hamiltonian matrix' diagonal. For the matrix elements of the kinetic energy operator, a band-limited (discrete $x_i$!) plane wave basis is used and the matrix elements $\left< x_i\vert T\vert x_j\right>$ are computed analytically.

You'll find a great discussion of this and various related Fourier methods in David Tannor's book "Introduction to Quantum Mechanics: a time-dependent perspective". Original references are also given there.

You'll find a Jupyter notebook with a simple Python/Cython implementation that tackles a quartic potential here. To change the potential you work with, find the line with pot = lambda x: x**4 - 20*x**2 and change it to any other potential you're interested in. Make sure to compute only bound states for which $\left<\psi\vert\psi\right>$ has decayed to zero close to the boundary of your $x$ domain.

The standard way to find the eigenvalues of the Schrodinger equation is called "imaginary time propagation". You change the coordinates, t=-i\tau, and integrate in the \tau direction. Any random initial condition will converge to the lowest energy eigenstate. The resulting equation is solved by splitting methods: First propagate the kinetic energy using Fast Fourier Transforms and then propagate the potential. This is possible since the solution decays fast (exponentially) as |x| grows and you can effectively replace your boundary conditions by periodic ones.

Here is some matlab code for a second order method

% in order to set up the problem, you need a spatial grid and a momentum grid,
% since you will be using FFTs, the number of grid points should be a power of 2
Dx  = 0.01;
x   = [-1:Dx:1-Dx];

% Now you define the eigenvalues in Fourier space
k           = zeros(1,Len).';
k(2:Len/2)  = 1:Len/2-1;                %  (x, 1, ..., Len/2-1, x, x, ..., x)
k((Len/2+1):Len) = -k((Len/2+1):-1:2);  %  (x, ...,x, -Len/2+1, ..., -1)
k(Len/2+1) = 0;

% from this you derive the momentum operator in Fourier space
p = -1i k.^2;

% This allows you to define the full problem, kinetic energy and potential, adjust as you please
E_kin = p.^2 / 2/mass;
V     = x.^4 - x.^2 ;

% you want to solve e^(h(T+V))Psi
% you approximate by the Strang splitting e^(h/2 V)e^(h T)e^(h/2 V)
for j=1:1:n
Psi = ifft( Psi);
Psi =  exp(- 1i*h* E_Kin ) .* Psi;
Psi =  fft( Psi);
% Apply potential V
Psi = exp(  -1i*h/2* V)
% Normalize the result in each step
Psi = Psi / ( norm(Psi,2)^2*Dx);
end
end


More sophisticated (higher order) methods can be found here, as well as many references to the problem: http://personales.upv.es/serblaza/2013JChemPhys.pdf

Since you want more eigenstates, there are several procedures: 1: iteration. compute the first eigenstate. Then, you restart the procedure, but in every step, you subtract the projection to this state. E.g. for the second state psi_2, in each step you write psi_2 = psi_2 - *psi_1.

As a complete alternative, you can discretize the Laplacian using finite differences (or use again FFTs) and then use the inverse power method. given your Hamiltonian H, iterate over n: (H-\lambda^(-1)I)^(-1)\psi_(n+1) = \psi_n for some guess of the eigenvalue \lambda. You can find details on wikipedia.

• Are you one of the author in the reference in the hyperlink? If so, you should explicitly state it. – nicoguaro Aug 9 '16 at 18:29
• why would that be necessary? – phil Aug 10 '16 at 2:02
• Is to avoid self promotion. That way people know that you are the author. – nicoguaro Aug 10 '16 at 2:21
• While imaginary time propagation is certainly an appealing method, I wouldn't call it the "standard way" of solving Schrödinger's equation. You may check e.g. Marston's Fourier Grid Hamiltonian Method (scitation.aip.org/content/aip/journal/jcp/91/6/10.1063/1.456888 or www.chem.yale.edu/~batista/classes/v572/fourier_grid_hamiltonian.pdf) for a simple yet effective&efficient method. – AlexE Aug 10 '16 at 18:38
• Maybe I should have specified that it is the standard method if the first few eigenvalue/eigenvector pairs are sought. The method referenced by AlexE solves the entire fully discretised Hamiltonian and makes sense only if many of such pairs are needed which didn't seem to be the case here (or possibly ever since there are infinitely many anyway). – phil Aug 11 '16 at 0:42