pyalps: Worm algorithm calculation of the two-component Bose-Hubbard mode

[witchbunting-cmd]

Hi！
I am attempting to use the worm algorithm in the pyalps library to calculate the two-component Bose-Hubbard model. The two-component Bose-Hubbard model implies that each lattice site can host two different types of atoms, and its Hamiltonian is as follows:

To calculate this model, I have constructed a model.xml file of the following form:

<SITEBASIS name="twoboson_site">
    <PARAMETER name="Na_max" default="2"/>
    <PARAMETER name="Nb_max" default="1"/>
    <QUANTUMNUMBER name="Na" min="0" max="Na_max"/>
    <QUANTUMNUMBER name="Nb" min="0" max="Nb_max"/>
    <!-- A -->
    <OPERATOR name="bdag_a" matrixelement="sqrt(Na+1)">
      <CHANGE quantumnumber="Na" change="1"/>
    </OPERATOR>
    <OPERATOR name="b_a" matrixelement="sqrt(Na)">
      <CHANGE quantumnumber="Na" change="-1"/>
    </OPERATOR>
    <OPERATOR name="nA" matrixelement="Na"/>
    <!-- B -->
    <OPERATOR name="bdag_b" matrixelement="sqrt(Nb+1)">
      <CHANGE quantumnumber="Nb" change="1"/>
    </OPERATOR>
    <OPERATOR name="b_b" matrixelement="sqrt(Nb)">
      <CHANGE quantumnumber="Nb" change="-1"/>
    </OPERATOR>
    <OPERATOR name="nB" matrixelement="Nb"/>
  </SITEBASIS>
  <BASIS name="twoboson_basis">
    <SITEBASIS ref="twoboson_site"/>
  </BASIS>
  <HAMILTONIAN name="twoboson">
    <PARAMETER name="J_a" default="1"/>
    <PARAMETER name="J_b" default="1"/>
    <PARAMETER name="U_a" default="0"/>
    <PARAMETER name="U_b" default="0"/>
    <PARAMETER name="U_ab" default="0"/>
    <PARAMETER name="mu_a" default="0"/>
    <PARAMETER name="mu_b" default="0"/>
    <BASIS ref="twoboson_basis"/>
    <SITETERM site="i">
      (U_a/2)*nA(i)*(nA(i)-1)
      + (U_b/2)*nB(i)*(nB(i)-1)
      + U_ab*nA(i)*nB(i)
      - mu_a*nA(i) - mu_b*nB(i)
    </SITETERM>
    <BONDTERM source="i" target="j">
      -J_a * (bdag_a(i)*b_a(j) + bdag_a(j)*b_a(i))
      -J_b * (bdag_b(i)*b_b(j) + bdag_b(j)*b_b(i))
    </BONDTERM>
  </HAMILTONIAN>

Then, the parameters are input and executed using the following Python code:

parms = []

parms.append({
        'LATTICE'          : "square lattice",    
        'MODEL'            : "twoboson",      
        'T'                : 0.125,                
        'J_a'              : 1.0,
        'J_b'              : 1.0,
        'U_a'              : 4.0,
        'U_b'              : 4.0,
        'U_ab'             : 2.0,
        'mu_a'             : 1.0,
        'mu_b'             : 1.0,           
        'Namax'            : 4,                 
        'Nbmax'            : 1,
        'L'                : 6,          

        'ALGORITHM'       : 'worm',             
        'THERMALIZATION'  : 10000,         
        'SWEEPS'          : 40000,          
        'MEASURE_AVERAGE[A_density]':"nA",   
        'MEASURE_AVERAGE[B_density]':"nB"        
    })
input_file = pyalps.writeInputFiles(dirs+'bh_worm', parms)
res = pyalps.runApplication('worm', input_file, writexml=True)
data = pyalps.loadMeasurements(pyalps.getResultFiles(dirs,prefix='bh_worm'), 
                               ['A_density', 'B_density','Energy'])

The code ran successfully and produced measurement results, but the average particle density nA is always zero, while nB outputs normally. Currently, it is unclear whether the issue lies with the measurement setup or if the algorithm does not support the two-component Hamiltonian on the same lattice site. If there are any errors, please point them out, and it would be greatly appreciated if you could provide suggestions for correction.

Thanks!

[LodePollet]

you could model your two-component system as a one-component system with two sheets of your two-dimensional square lattice where the "A" particles live on the lower sheet and the "B" particles on the upper sheet.