1function [OP,P1m]= QuantumSimulationCavityArray(wc,D,g,Nphoton,L,A,Sp,N_ex,Hhopp,t,Model)
2HJC = zeros((Nphoton+1)^L*2^L,(Nphoton+1)^L*2^L);
3for i=1:L
4 switch Model
5 case 'Jaynes-Cummings'
6 HJC = HJC + wc*A{i}'*A{i} + (D+wc)*Sp{i}*Sp{i}' + g*(Sp{i}*A{i}+Sp{i}'*A{i}');
7 case 'Rabi'
8 HJC = HJC + wc*A{i}'*A{i} + (D+wc)*Sp{i}*Sp{i}' + g*(Sp{i}+Sp{i}')*(A{i}+A{i}');
9 end
10end
11H = HJC + Hhopp; % Total Hamiltonian
12dt = t(2)-t(1); % Step time dt
13U = expm(-1i*H*dt); % Time propagator with step dt
14up = [1 0]'; % Excited state for the two-level system
15down = [0 1]'; % Ground state for the two-level system
16Fock = eye(Nphoton+1); % Fock states
17theta1 = 0.5*atan(2*g*sqrt(1)/D);
18phi_1m = cos(theta1)*kron(down,Fock(:,2))...
19 -sin(theta1)*kron(up,Fock(:,1)); % Initial state |1,->
20PSI_0 = phi_1m;
21for k=1:L-1
22 PSI_0 = kron(PSI_0,phi_1m); % Many body wavefunction |1,->...|1,-> (L terms)
23end
24dn_T = zeros(size(t)); % Standard deviation dn_i = <n_i^2>-<n_i>^2
25P1m = zeros(size(t)); % Population |1,-> at time t
26for n=1:length(t)
27 if n==1
28 PSI = PSI_0; % Initial wavefunction
29 else
30 PSI = U*PSI; % Wavefunction at time t_n
31 end
32 dn = 0;
33 for i=1:L
34 dn = dn + PSI'*N_ex{i}^2*PSI-(PSI'*N_ex{i}*PSI)^2;
35 end
36 P1m(n) = abs(PSI_0'*PSI)^2; % Ground state probability
37 dn_T(n) = dn; % Standard deviation of the number of excitations
38end
39OP = mean(dn_T); % Order parameter
40end