How to split objects to increase effective resolution of particle effects
As the particle system is just a way to obtain the result, I propose to do it another way. A script was needed anyway, so the following script does a voxelisation directly.



Instead of PS, it uses a script and dupliverts.

The script creates a mesh with vertices only so that the duplivert can be applied:


It takes a mesh as input
Calculates its boundings
Maps x, y, z inside the boundings at voxel size steps
Then uses a BVHTree to test the distance of these (x, y, z) to the original mesh


Once done, you just have to create a cube at the same voxel size and parent it to the generated mesh (with duplivert).

From that, you may want to make an exportable mesh:


Select the dupli
Use the menu 'object>apply>make instance real'
Select the instances and join them


Performance: tested on this bunny at 2mm voxel size. That's about 350k vertices. So 1mm should be around 8 times this amount.

The calculation is relatively fast, but once the duplivert set up, Blender becomes a bit slow.

The script:

import bpy
import bmesh
from math import ceil
from mathutils import Vector
from mathutils.bvhtree import BVHTree
from itertools import product


# Create an object from vertices only
def create_object( context, vertices ):
    mesh = bpy.data.meshes.new( "voxel" )
    obj = bpy.data.objects.new( "voxel", mesh )
    context.scene.collection.objects.link( obj )
    bm = bmesh.new()

    for v in vertices:
        bm.verts.new(v)

    bm.to_mesh( mesh )
    bm.free()

    return obj

# Get object bounding as min max vectors
def boundings( obj ):
    min_x = min( (Vector(b).x for b in obj.bound_box) )
    max_x = max( (Vector(b).x for b in obj.bound_box) )
    min_y = min( (Vector(b).y for b in obj.bound_box) )
    max_y = max( (Vector(b).y for b in obj.bound_box) )
    min_z = min( (Vector(b).z for b in obj.bound_box) )
    max_z = max( (Vector(b).z for b in obj.bound_box) )
    return Vector( (min_x, min_y, min_z) ), Vector( (max_x, max_y, max_z) ) 

# Voxel
def voxelize( context, obj, size ):
    # Get boundings
    min_bounds, max_bounds = boundings( obj )
    # Delta vector
    delta_bounds = max_bounds - min_bounds
    # Vertex count in all dimensions
    counts = Vector( (ceil( delta_bounds.x / size ), ceil( delta_bounds.y / size ), ceil( delta_bounds.z / size )) )
    # Reajusted starting position
    starts = min_bounds - ((counts * size) - delta_bounds) * 0.5

    # A BVHTree from the mesh
    depsgraph = context.evaluated_depsgraph_get()
    tree = BVHTree.FromObject( obj, depsgraph )

    half_size = size / 2.0
    vertices = [] 
    # Cross products of the wanted ranges
    for x, y, z in product( range(int(counts.x)+1), range(int(counts.y)+1), range(int(counts.z)+1) ):
        # Corresponding position
        vert = starts + Vector( (x * size, y * size, z * size) )

        # Test the proximity to the mesh with the BVH tree
        result = tree.find_nearest_range( vert, half_size )
        # If close enough add the location
        if result:
            vertices.append( vert )

    # Create an mesh from the result
    return create_object( context, vertices )


# Get current object    
obj = bpy.context.object
# Define voxel size
voxel_size = 0.002 # &lt;&lt;&lt; Change this to change the resolution (in meters)

# Create the voxel
voxelize( bpy.context, obj, voxel_size )


The file with low res voxel:



Another example:



Complements:

To use the script:


This is not an addon to Blender, so you have to open a text view in Blender, create a new text and paste the code in it (same way it is in the uploaded blend file)
The script contains a parameter in the line commented '&lt;&lt;&lt; Change this to change the resolution (in meters)' above. This 'voxel_size' parameter corresponds to the size of the cubes
Then select the mesh (active mesh) in the 3D View and run the script using the button 'run script' which is on the top right of the menus of the text view (you can also use AltP when the mouse is in the text view)
The calculation starts after that and Blender may hang for some seconds. This calculation time depends on the 'voxel_size' parameter and on the mesh dimensions. The ratio 'dimensions/voxel_size' power 3 gives the voxel cubes amount. I preconize you start with a relatively high value for voxel_size so that there is few cubes in result. After that and depending on your configuration lower this value.
Once done, the script has create a new mesh which is vertex only (no faces). You'll have to follow what is written above about dupli and parenting to make a base cube of voxel_size dimension fill the mesh