I think it's best to do this in two steps. In the first step, we rotate the bone to align its y-axis with your custom vector. Then, in the second step, we rotate the bone along its y-axis (or the custom vector) to align its x-axis with its parent's x-axis. To do this, we'll need to project the parent's x-axis to the xz-plane of the bone. The following screenshots illustrate the two steps.
I decided to do this in the local space of the armature, so I needed to add an additional step to account for the armature object rotation. If the armature object is aligned with the world, the extra step can be skipped.
The script is attached below with comments explaining how I did it.
import bpy, mathutils
# Define the target vector
v = mathutils.Vector((0.0, 0.0, 1.0))
# Get the armature object by object name
arm = bpy.data.objects['Armature']
# Get the active bone. If you want to operator on a specific bone, use it's name instead. arm.pose.bones['Bone']
bone = arm.pose.bones[arm.data.bones.active.name]
# Get the components of the matrix of the bone. matrix is in the armature object space. matrix_basis is relative to the parent. decompose() returns the location, rotation, and scale from the matrix.
l, r, s = bone.matrix.decompose()
# The first step:
# rotation_difference returns a quaternion representing the rotation between two vectors. In our case, we need a quaternion that represent a rotation to our target vector.
r1 = bone.y_axis.rotation_difference(v)
# Apply the rotation to the bone's rotation
r.rotate(r1)
# The extra step:
# To cancel out the armature object rotation, we need to extract the rotation from its matrix_world and apply its inverse to the bone's rotation.
r_extra = arm.matrix_world.decompose()[1]
r.rotate(r_extra.inverted())
# The second step depends on our work so far being applied. We do that as follows.
# Update the matrix of the bone by building and assigning a new matrix from the old location and scale, and the modified rotation
bone.matrix = mathutils.Matrix.LocRotScale(l, r, s)
# Trigger an update.
bpy.context.view_layer.update()
# The second step:
# Repeating the first step with the parent's x-axis and the bone's x-axis will rotate the whole bone. This is not what we want. We want the bone to rotate about its y-axis. This means we need to use the parent's x-axis projection on the bone's xz-plane instead.
# I couldn't find a direct way to achieve this in the API, so I constructed the projection with the next best thing I could find.
parent_x = bone.parent.x_axis
x_proj = parent_x.project(bone.x_axis)
z_proj = parent_x.project(bone.z_axis)
parent_x_proj = x_proj + z_proj
# The remaining steps were explained previously
r2 = bone.x_axis.rotation_difference(parent_x_proj)
r.rotate(r2)
bone.matrix = mathutils.Matrix.LocRotScale(l, r, s)
